text:
using Molecules.PointGroup

mols = ["water", "ammonia", "benzene", "ethanol", "formaldehyde", "glycine", "methane", "phosphaethene", "acetic_acid", "nitrogen"]
Ds = Any[Molecules.PointGroup.DistanceMatrix(Any[8, 1, 1], [0.0 0.9902542038608343 0.9902542038706021; 0.9902542038608343 0.0 1.5660733799622173; 0.9902542038706021 1.5660733799622173 0.0]), Molecules.PointGroup.DistanceMatrix(Any[7, 1, 1, 1], [0.0 1.0113541286498018 1.011354577368887 1.011354577368887; 1.0113541286498018 0.0 1.6153960232339313 1.6153960232339313; 1.011354577368887 1.6153960232339313 0.0 1.615396; 1.011354577368887 1.6153960232339313 1.615396 0.0]), Molecules.PointGroup.DistanceMatrix(Any[6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1], [0.0 1.391673101250793 2.4104478246286516 2.783346 2.4104478246286516 1.391673101250793 1.077915 2.1444719451923357 3.3871046345167724 3.861261 3.3871046345167724 2.1444719451923357; 1.391673101250793 0.0 1.391672 2.4104478246286516 2.7833452025014793 2.410448 2.144472207252871 1.0779154464836285 2.1444715388421454 3.387104071319923 3.861260648984992 3.3871045641172643; 2.4104478246286516 1.391672 0.0 1.391673101250793 2.410448 2.7833452025014793 3.387104071319923 2.1444715388421454 1.0779154464836285 2.144472207252871 3.3871045641172643 3.861260648984992; 2.783346 2.4104478246286516 1.391673101250793 0.0 1.391673101250793 2.4104478246286516 3.861261 3.3871046345167724 2.1444719451923357 1.077915 2.1444719451923357 3.3871046345167724; 2.4104478246286516 2.7833452025014793 2.410448 1.391673101250793 0.0 1.391672 3.387104071319923 3.861260648984992 3.3871045641172643 2.144472207252871 1.0779154464836285 2.1444715388421454; 1.391673101250793 2.410448 2.7833452025014793 2.4104478246286516 1.391672 0.0 2.144472207252871 3.3871045641172643 3.861260648984992 3.387104071319923 2.1444715388421454 1.0779154464836285; 1.077915 2.144472207252871 3.387104071319923 3.861261 3.387104071319923 2.144472207252871 0.0 2.4695880477342778 4.277451917321806 4.939176 4.277451917321806 2.4695880477342778; 2.1444719451923357 1.0779154464836285 2.1444715388421454 3.3871046345167724 3.861260648984992 3.3871045641172643 2.4695880477342778 0.0 2.469588 4.277451917321806 4.9391760954685555 4.277452; 3.3871046345167724 2.1444715388421454 1.0779154464836285 2.1444719451923357 3.3871045641172643 3.861260648984992 4.277451917321806 2.469588 0.0 2.4695880477342778 4.277452 4.9391760954685555; 3.861261 3.387104071319923 2.144472207252871 1.077915 2.144472207252871 3.387104071319923 4.939176 4.277451917321806 2.4695880477342778 0.0 2.4695880477342778 4.277451917321806; 3.3871046345167724 3.861260648984992 3.3871045641172643 2.1444719451923357 1.0779154464836285 2.1444715388421454 4.277451917321806 4.9391760954685555 4.277452 2.4695880477342778 0.0 2.469588; 2.1444719451923357 3.3871045641172643 3.861260648984992 3.3871046345167724 2.1444715388421454 1.0779154464836285 2.4695880477342778 4.277452 4.9391760954685555 4.277451917321806 2.469588 0.0]), Molecules.PointGroup.DistanceMatrix(Any[6, 6, 8, 1, 1, 1, 1, 1, 1], [0.0 1.5066066187356275 2.3566868414861575 3.1934221899825586 1.0864675442708818 1.08568838951515 1.08568838951515 2.135170299478709 2.135170299478709; 1.5066066187356275 0.0 1.4125683427997386 1.9398012042879549 2.143783442067785 2.1382857296381137 2.1382857296381137 1.0907794211709352 1.0907794211709352; 2.3566868414861575 1.4125683427997386 0.0 0.9547383645423493 3.3084015189628966 2.606508193625717 2.606508193625717 2.073565336203805 2.073565336203805; 3.1934221899825586 1.9398012042879549 0.9547383645423493 0.0 4.042906164996289 3.4835114492754293 3.4835114492754293 2.3207982562198293 2.3207982562198293; 1.0864675442708818 2.143783442067785 3.3084015189628966 4.042906164996289 0.0 1.7663523053700245 1.7663523053700245 2.466529679033682 2.466529679033682; 1.08568838951515 2.1382857296381137 2.606508193625717 3.4835114492754293 1.7663523053700245 0.0 1.762268 2.4748288131844594 3.0380137127445623; 1.08568838951515 2.1382857296381137 2.606508193625717 3.4835114492754293 1.7663523053700245 1.762268 0.0 3.0380137127445623 2.4748288131844594; 2.135170299478709 1.0907794211709352 2.073565336203805 2.3207982562198293 2.466529679033682 2.4748288131844594 3.0380137127445623 0.0 1.761792; 2.135170299478709 1.0907794211709352 2.073565336203805 2.3207982562198293 2.466529679033682 3.0380137127445623 2.4748288131844594 1.761792 0.0]), Molecules.PointGroup.DistanceMatrix(Any[8, 6, 1, 1], [0.0 1.2042169999999999 2.0134166958799167 2.0134166958799167; 1.2042169999999999 0.0 1.0995655003723062 1.0995655003723062; 2.0134166958799167 1.0995655003723062 0.0 1.869456; 2.0134166958799167 1.0995655003723062 1.869456 0.0]), Molecules.PointGroup.DistanceMatrix(Any[6, 8, 8, 6, 7, 1, 1, 1, 1, 1], [0.0 1.2062396067999923 1.351799423276989 1.5236746177120626 2.5105614991847935 1.8654645741809195 2.1370666272877408 2.1370666272877408 2.668320812784887 2.668320812784887; 1.2062396067999923 0.0 2.2524002769705476 2.4267294659019987 2.847682130443986 2.2784360771039855 3.1205223347237876 3.1205223347237876 2.7193182166969723 2.7193182166969723; 1.351799423276989 2.2524002769705476 0.0 2.380395481079562 3.6846896878468614 0.9703246445350132 2.607391363669635 2.607391363669635 3.9135222900452478 3.9135222900452478; 1.5236746177120626 2.4267294659019987 2.380395481079562 0.0 1.4556534497334865 3.2031556140843356 1.1051975920906632 1.1051975920906632 2.017632237009758 2.017632237009758; 2.5105614991847935 2.847682130443986 3.6846896878468614 1.4556534497334865 0.0 4.368948064256772 2.1099113358788326 2.1099113358788326 1.02388218454957 1.02388218454957; 1.8654645741809195 2.2784360771039855 0.9703246445350132 3.2031556140843356 4.368948064256772 0.0 3.4865087194319475 3.4865087194319475 4.472711679825673 4.472711679825673; 2.1370666272877408 3.1205223347237876 2.607391363669635 1.1051975920906632 2.1099113358788326 3.4865087194319475 0.0 1.76263 2.4146608549943407 2.9431608409127086; 2.1370666272877408 3.1205223347237876 2.607391363669635 1.1051975920906632 2.1099113358788326 3.4865087194319475 1.76263 0.0 2.9431608409127086 2.4146608549943407; 2.668320812784887 2.7193182166969723 3.9135222900452478 2.017632237009758 1.02388218454957 4.472711679825673 2.4146608549943407 2.9431608409127086 0.0 1.606468; 2.668320812784887 2.7193182166969723 3.9135222900452478 2.017632237009758 1.02388218454957 4.472711679825673 2.9431608409127086 2.4146608549943407 1.606468 0.0]), Molecules.PointGroup.DistanceMatrix(Any[6, 1, 1, 1, 1], [0.0 1.085807062807661 1.085807062807661 1.085807062807661 1.085807062807661; 1.085807062807661 0.0 1.7731155086592638 1.7731155086592638 1.7731155086592638; 1.085807062807661 1.7731155086592638 0.0 1.7731155086592638 1.7731155086592638; 1.085807062807661 1.7731155086592638 1.7731155086592638 0.0 1.7731155086592638; 1.085807062807661 1.7731155086592638 1.7731155086592638 1.7731155086592638 0.0]), Molecules.PointGroup.DistanceMatrix(Any[6, 15, 1, 1, 1], [0.0 1.681005 1.0829540195968617 1.0846240077672076 2.335469999487255; 1.681005 0.0 2.462851028563441 2.4055647894517413 1.4254209965129598; 1.0829540195968617 2.462851028563441 0.0 1.839285554824481 2.528008098146048; 1.0846240077672076 2.4055647894517413 1.839285554824481 0.0 3.35847532580634; 2.335469999487255 1.4254209965129598 2.528008098146048 3.35847532580634 0.0]), Molecules.PointGroup.DistanceMatrix(Any[6, 6, 8, 1, 1, 1, 8, 1], [0.0 1.4919159841857719 2.390674944972863 1.078682443370615 1.0861889075943465 1.0861889075943465 2.360210142020833 3.1526672427076727; 1.4919159841857719 0.0 1.20495230257882 2.0839124584744435 2.1304824072763426 2.1304824072763426 1.3373790462053006 1.8240597581704938; 2.390674944972863 1.20495230257882 0.0 2.4945753654457903 3.0929792987360907 3.0929792987360907 2.2301024588599065 2.2171180459057656; 1.078682443370615 2.0839124584744435 2.4945753654457903 0.0 1.7772100119594758 1.7772100119594758 3.249353594407663 3.901379292748912; 1.0861889075943465 2.1304824072763426 3.0929792987360907 1.7772100119594758 0.0 1.757266 2.6486817262908358 3.5014817176508286; 1.0861889075943465 2.1304824072763426 3.0929792987360907 1.7772100119594758 1.757266 0.0 2.6486817262908358 3.5014817176508286; 2.360210142020833 1.3373790462053006 2.2301024588599065 3.249353594407663 2.6486817262908358 2.6486817262908358 0.0 0.9605159079786239; 3.1526672427076727 1.8240597581704938 2.2171180459057656 3.901379292748912 3.5014817176508286 3.5014817176508286 0.9605159079786239 0.0]), Molecules.PointGroup.DistanceMatrix(Any[7, 7], [0.0 1.100592; 1.100592 0.0])]
SEAs = Any[[[1], [2, 3]], [[1], [2, 3, 4]], [[1, 2, 3, 4, 5, 6], [7, 8, 9, 10, 11, 12]], [[1], [2], [3], [4], [5], [6, 7], [8, 9]], [[1], [2], [3, 4]], [[1], [2], [3], [4], [5], [6], [7, 8], [9, 10]], [[1], [2, 3, 4, 5]], [[1], [2], [3, 5], [4, 5]], [[1], [2], [3], [4], [5, 6], [7], [8]], [[1, 2]]]
COMs = Any[[1.2483188267782848, 1.8068607904101412, 0.02067611111507077], [0.0, -5.9186189890467085e-8, 0.04789390012330456], [0.0, -5.685142354124265e-18, 0.0], [-0.0541194764592242, -0.017770381536391034, 0.0], [0.0, 0.0, 0.07304868310797315], [0.012654881758962023, 0.040814231526502945, 0.0], [0.0, 0.0, 0.0], [0.02686983463465735, -0.06696346354690948, 0.0], [-0.030964198561246933, 0.05292127193099316, 0.0], [0.0, 0.0, 0.0]]
COMshifts = Any[Any[Atom{Int16, Float64}(8, 15.999, [-0.03916517197828484, -0.04044897151014126, -0.03783750831507077]), Atom{Int16, Float64}(1, 1.008, [0.9501611807217154, -0.009150727710141249, -0.00855993921507077]), Atom{Int16, Float64}(1, 1.008, [-0.3285306385782848, 0.6511577665898587, 0.6091177720849292])], Any[Atom{Int16, Float64}(7, 14.007, [0.0, 5.9186189890467085e-8, 0.06945309987669546]), Atom{Int16, Float64}(1, 1.008, [0.0, 0.9326490591861898, -0.3217029001233046]), Atom{Int16, Float64}(1, 1.008, [0.807698, -0.4663249408138101, -0.3217029001233046]), Atom{Int16, Float64}(1, 1.008, [-0.807698, -0.4663249408138101, -0.3217029001233046])], Any[Atom{Int16, Float64}(6, 12.011, [0.0, 1.391673, 0.0]), Atom{Int16, Float64}(6, 12.011, [1.205224, 0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [1.205224, -0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.0, -1.391673, 0.0]), Atom{Int16, Float64}(6, 12.011, [-1.205224, -0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [-1.205224, 0.695836, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.0, 2.469588, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.138726, 1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.138726, -1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.0, -2.469588, 0.0]), Atom{Int16, Float64}(1, 1.008, [-2.138726, -1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [-2.138726, 1.234794, 0.0])], Any[Atom{Int16, Float64}(6, 12.011, [1.2157024764592241, -0.38898461846360893, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.0541194764592242, 0.570488381536391, 0.0]), Atom{Int16, Float64}(8, 15.999, [-1.132994523540776, -0.19508961846360895, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.8783145235407759, 0.40158738153639106, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.1569794764592243, 0.15361038153639103, 0.0]), Atom{Int16, Float64}(1, 1.008, [1.176466476459224, -1.022058618463609, 0.881134]), Atom{Int16, Float64}(1, 1.008, [1.176466476459224, -1.022058618463609, -0.881134]), Atom{Int16, Float64}(1, 1.008, [0.1102664764592242, 1.211323381536391, 0.880896]), Atom{Int16, Float64}(1, 1.008, [0.1102664764592242, 1.211323381536391, -0.880896])], Any[Atom{Int16, Float64}(8, 15.999, [0.0, 0.0, 0.6014443168920268]), Atom{Int16, Float64}(6, 12.011, [0.0, 0.0, -0.6027726831079732]), Atom{Int16, Float64}(1, 1.008, [0.0, 0.934728, -1.1818476831079732]), Atom{Int16, Float64}(1, 1.008, [0.0, -0.934728, -1.1818476831079732])], Any[Atom{Int16, Float64}(6, 12.011, [-0.012654881758962023, 0.509800768473497, 0.0]), Atom{Int16, Float64}(8, 15.999, [1.164952118241038, 0.771058768473497, 0.0]), Atom{Int16, Float64}(8, 15.999, [-0.982078881758962, 1.451914768473497, 0.0]), Atom{Int16, Float64}(6, 12.011, [-0.597416881758962, -0.8971952315265029, 0.0]), Atom{Int16, Float64}(7, 14.007, [0.388003118241038, -1.968584231526503, 0.0]), Atom{Int16, Float64}(1, 1.008, [-0.519754881758962, 2.305018768473497, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.258313881758962, -0.986418231526503, 0.881315]), Atom{Int16, Float64}(1, 1.008, [-1.258313881758962, -0.986418231526503, -0.881315]), Atom{Int16, Float64}(1, 1.008, [1.005802118241038, -1.822043231526503, 0.803234]), Atom{Int16, Float64}(1, 1.008, [1.005802118241038, -1.822043231526503, -0.803234])], Any[Atom{Int16, Float64}(6, 12.011, [0.0, 0.0, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.626891, 0.626891, 0.626891]), Atom{Int16, Float64}(1, 1.008, [-0.626891, -0.626891, 0.626891]), Atom{Int16, Float64}(1, 1.008, [-0.626891, 0.626891, -0.626891]), Atom{Int16, Float64}(1, 1.008, [0.626891, -0.626891, -0.626891])], Any[Atom{Int16, Float64}(6, 12.011, [0.02979516536534265, 1.1472954635469095, 0.0]), Atom{Int16, Float64}(15, 30.9737619985, [0.02979516536534265, -0.5337095364530905, 0.0]), Atom{Int16, Float64}(1, 1.008, [-0.8617048346346573, 1.7621264635469094, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.9756601653653426, 1.6780954635469094, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.3845288346346574, -0.7112275364530904, 0.0])], Any[Atom{Int16, Float64}(6, 12.011, [1.0153271985612469, -1.0325762719309932, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.030964198561246933, 0.08851372806900684, 0.0]), Atom{Int16, Float64}(8, 15.999, [0.32640719856124695, 1.2566847280690068, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.000327198561247, -0.5928822719309932, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.8801711985612469, -1.6567102719309932, 0.878633]), Atom{Int16, Float64}(1, 1.008, [0.8801711985612469, -1.6567102719309932, -0.878633]), Atom{Int16, Float64}(8, 15.999, [-1.238545801438753, -0.33211327193099316, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.7504518014387531, 0.4806247280690068, 0.0])], Any[Atom{Int16, Float64}(7, 14.007, [0.0, 0.0, 0.550296]), Atom{Int16, Float64}(7, 14.007, [0.0, 0.0, -0.550296])]]
MOITs = Any[[59.67273019851896 -40.43446779974129 -0.2787699419816882; -40.43446779974129 29.520811530245208 -1.0973898240739823; -0.2787699419816882 -1.0973898240739823 88.38419501600158], [1.734784454040015 0.0 0.0; 0.0 1.734784403596671 -2.759994720546999e-7; 0.0 -2.759994720546999e-7 2.630380321153872], [88.2299814706577 0.0 0.0; 0.0 88.22999596081817 0.0; 0.0 0.0 176.45997743147586], [14.729848140524304 4.309547442732232 0.0; 4.309547442732232 52.64956324357207 0.0; 0.0 0.0 61.12023583256766], [14.888928588039047 0.0 0.0; 0.0 13.127516257127304 0.0; 0.0 0.0 1.761412330911744], [127.31198181986389 15.141326500943165 0.0; 15.141326500943165 51.92297634347139 0.0; 0.0 0.0 173.50185337979588], [3.1690901159043836 0.0 0.0; 0.0 3.1690901159043836 0.0; 0.0 0.0 3.1690901159043836], [31.31733220246329 -0.9476165389488573 0.0; -0.9476165389488573 3.7116365967505667 0.0; 0.0 0.0 35.028968799213864], [47.77658951957811 4.497413524493793 0.0; 4.497413524493793 48.93811681454731 0.0; 0.0 0.0 93.60201866901137], [8.483358812874624 0.0 0.0; 0.0 8.483358812874624 0.0; 0.0 0.0 0.0]]
EMOITs = Any[([1.4302678930357826, 87.5617576664536, 88.58571118527635], [0.5702815858999435 0.7449168596591891 0.346233714388038; 0.8213587619322918 -0.5108359683198183 -0.2538038566842631; 0.012194137136251365 -0.42912176088640136 0.9031643357403217]), ([1.7347844035965863, 1.734784454040015, 2.6303803211539574], [0.0 1.0 0.0; 0.9999999999999525 0.0 -3.081741069203038e-7; 3.081741069203038e-7 0.0 0.9999999999999525]), ([88.2299814706577, 88.22999596081817, 176.45997743147586], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]), ([14.24623898635917, 53.1331723977372, 61.12023583256766], [0.9937624015291687 0.11151811201315731 0.0; -0.11151811201315731 0.9937624015291688 0.0; 0.0 0.0 1.0]), ([1.761412330911744, 13.127516257127304, 14.888928588039047], [0.0 0.0 1.0; 0.0 1.0 0.0; 1.0 0.0 0.0]), ([48.99562160458952, 130.23933655874578, 173.50185337979588], [-0.1898203588075834 0.9818188383720086 0.0; 0.9818188383720086 0.18982035880758344 0.0; 0.0 0.0 1.0]), ([3.1690901159043836, 3.1690901159043836, 3.1690901159043836], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]), ([3.679146146650979, 31.349822652562874, 35.028968799213864], [0.03426636016583991 0.9994127358408962 0.0; 0.9994127358408962 -0.03426636016583991 0.0; 0.0 0.0 1.0]), ([43.82259685337404, 52.89210948075138, 93.60201866901137], [-0.7510224457563778 0.660276673804329 0.0; 0.6602766738043291 0.751022445756378 0.0; 0.0 0.0 1.0]), ([0.0, 8.483358812874624, 8.483358812874624], [0.0 0.0 1.0; 0.0 1.0 0.0; 1.0 0.0 0.0])]
tol = 1E-5

@testset "SymTools" begin
    for i = eachindex(mols)
        path = joinpath(@__DIR__, "xyz/"*mols[i]*".xyz")
        mol = Molecules.parse_file(path)

        D = buildD(mol)
        @test isapprox(D.distances, Ds[i].distances, rtol=tol)

        sea = findSEA(D, 5)
        @test isapprox(sea, SEAs[i], rtol=tol)

        moit = calcmoit(mol)
        @test isapprox(moit, MOITs[i], rtol=tol)

        emoit = eigenmoit(moit)
        @test isapprox(emoit[1], EMOITs[i][1], rtol=tol)
        @test isapprox(emoit[2], EMOITs[i][2], rtol=tol)

    end
end


text:
md"""
# Quick Tutorial

## Construct a very simple graph
Just to get started, lets create a simple graph for the summation of two numbers.
"""
using NaiveNASlib, Test
md"""
NaiveNASlib uses a special immutable type of vertex to annotate inputs so that one
can be certain that size mutations won't suddenly change the input shape of the model.
"""
@testset "First example" begin #src
in1 = inputvertex("in1", 1)
in2 = inputvertex("in2", 1)

# In this example we could have done without them as we won't have any parameters to
# change the size of, but lets show things as they are expected to be used.
#
# Create a new vertex which computes the sum of `in1` and `in2`:
add = "add" >> in1 + in2
@test add isa NaiveNASlib.AbstractVertex

# NaiveNASlib lets you do this using [`+`](@ref) which creates a new vertex which sums it inputs.
# Use `>>` to attach a name to the vertex when using infix operations.
# Naming vertices is completely optional, but it is quite helpful as can be seen below.

# [`CompGraph`](@ref) helps evaluating the whole graph as a function.
graph = CompGraph([in1, in2], add)

# Evaluate the function represented by `graph` by just calling it.
@test graph(2,3) == 5
@test graph(100,200) == 300

# The [`vertices`](@ref) function returns the vertices in topological order.
@test vertices(graph) == [in1, in2, add]
@test name.(vertices(graph)) == ["in1", "in2", "add"]

# [`CompGraph`](@ref)s can be indexed:
@test graph[begin] == graph[1] == in1
@test graph[end] == graph[3] == add
@test graph[begin:end] == vertices(graph)
# This is a bit slow though as it traverses the whole graph each time. It is better to
# call [`vertices`](@ref) first and then apply the indexing if one needs to do this many times.
end #src

# ## Modify a graph
# Now lets look at how to make use of it to modify the structure of a neural network. Since batteries are excluded,
# lets first create a tiny neural network library and do the wiring so that NaiveNASlib can work with it.

module TinyNNlib
    using NaiveNASlib
    ## A simple linear layer
    mutable struct LinearLayer{T}
        W::Matrix{T}
    end
    ## Normally ones uses something like randn here, but this makes output
    ## in examples easier on the eyes
    LinearLayer(nin, nout) = LinearLayer(ones(Int, nout, nin))
    (l::LinearLayer)(x) = l.W * x

    ## NaiveNASlib needs to know what LinearLayer considers its output and input size
    ## In this case it is the number of rows and columns of the weight matrix
    ## Input size is always a vector since vertices might have multiple inputs
    NaiveNASlib.nin(l::LinearLayer) = [size(l.W, 2)]
    NaiveNASlib.nout(l::LinearLayer) = size(l.W, 1)

    ## We also need to tell NaiveNASlib how to change the size of LinearLayer
    ## The Δsize! function will receive indices to keep from existing weights
    ## as well as where to insert new indices
    function NaiveNASlib.Δsize!(l::LinearLayer, newins::AbstractVector, newouts::AbstractVector)
        ## newins is a vector of vectors as vertices may have more than one input,
        ## but LinearLayer has only one
        ## The function NaiveNASlib.parselect can be used to interpret newins and newouts.
        ## We just need to tell it along which dimensions to apply them.
        l.W = NaiveNASlib.parselect(l.W, 1=>newouts, 2=>newins[])
    end

    ## Helper function which creates a LinearLayer wrapped in an vertex in a computation graph.
    ## This creates a Keras-like API
    linearvertex(in, outsize) = absorbvertex(LinearLayer(nout(in), outsize), in)
    export linearvertex, LinearLayer
end

md"""
There are a handful of other functions one can implement to e.g. provide better defaults and offer other forms of convenience,
but here we use the bare minimum to get ourselves started. Some examples of this is provided further down.

In practice one might not want to create a whole neural network library from scratch, but rather incorporate NaiveNASlib with an existing library
in a glue package.
It might appear that this will inevitably lead to type-piracy as neither the layer definitions nor the functions (e.g. nout, nin) would
belong to the glue package. However, one typically wants to wrap the layers in some intermediate type anyways. For instance,
[NaiveNASflux](https://github.com/DrChainsaw/NaiveNASflux.jl) needs to wrap the layers from Flux in a mutable container as the layers
themselves are not mutable.

Lets do a super simple example where we make use of the tiny neural network library to create a model and then modify it:
"""

@testset "Second example" begin #src
using .TinyNNlib
invertex = inputvertex("input", 3)
layer1 = linearvertex(invertex, 4)
layer2 = linearvertex(layer1, 5)

# Vertices may be called to execute their computation alone.
# We generally outsource this work to [`CompGraph`](@ref), but now we are trying to illustrate how things work.
batchsize = 2
batch = randn(nout(invertex), batchsize)
y1 = layer1(batch)
@test size(y1) == (nout(layer1), batchsize) == (4, 2)
y2 = layer2(y1)
@test size(y2) == (nout(layer2), batchsize) == (5, 2)

# Lets change the output size of `layer1`. First check the input sizes so we have something to compare to.
@test [nout(layer1)] == nin(layer2) == [4]
@test Δnout!(layer1 => -2) # Returns true if successful
@test [nout(layer1)] == nin(layer2) == [2]

# And now the weight matrices have changed!
#
# The graph is still operational of course but the sizes of the activations have changed.
y1 = layer1(batch)
@test size(y1) == (nout(layer1), batchsize) == (2, 2)
y2 = layer2(y1)
@test size(y2) == (nout(layer2), batchsize) == (5 ,2)
end #src

md"""

As can be seen above, the consequence of changing the output size of `layer1` was that the input size of `layer2`
also was changed. This is of course required for the computation graph to not throw a dimension mismatch error
when being executed.

Besides the very simple graph, this mutation was trivial because the sizes of the input and output dimensions
of `LinearLayer`'s parameters can change independently as they are the rows and columns of the weight matrix.
This is expressed by giving the layers the mutation size trait [`SizeAbsorb`](@ref) in the lingo of NaiveNASlib,
meaning that a change in number of input/output neurons does not propagate further in the graph.

## A more elaborate example

While the previous example was simple enough to be done by hand, things can quickly get out of hand when using:

* Layers which require `nin==nout`, e.g. batch normalization and pooling.
* Element wise operations such as activation functions or just element wise arithmetics (e.g `+` used in residual connections).
* Concatenation of activations.

Lets use a small but non-trivial model including all of the above. We begin by making a helper which creates a vertex which
does elementwise scaling of its input:
"""
scalarmult(v, s::Number) = invariantvertex(x -> x .* s, v)
# When multiplying with a scalar, the output size is the same as the input size.
# This vertex type is said to be size invariant (in lack of better words), hence the name [`invariantvertex`](@ref).

@testset "More elaborate example" begin #src
# Ok, lets create the model:
## First a few `LinearLayer`s
invertex = inputvertex("input", 6)
start = linearvertex(invertex, 6)
split = linearvertex(start, nout(invertex) ÷ 3)

## Concatenation means the output size is the sum of the input sizes
joined = conc(scalarmult(split,2), scalarmult(split,3), scalarmult(split,5), dims=1)

## Elementwise addition is of course also size invariant
out = start + joined

## CompGraph to help us run the whole thing
graph = CompGraph(invertex, out)
@test graph((ones(6)))  == [78, 78, 114, 114, 186, 186]

md"""
Now we have a somewhat complex set of size relations at our hand since the sizes are constrained so that
  1. `start` and `joined` must have the same output size due to element wise addition.
  2. `joined` will always have 3 times the output size of `split` since there are no size absorbing vertices between them.

Modifying this graph manually would of course be manageable (albeit a bit cumbersome) if we created the model
by hand and knew it inside out. When things like the above emerges out of a neural architecture search things
are less fun though and this is where use of NaiveNASlib will really pay off.

Ok, lets try to increase the size of the vertex `out` by 2. Before we do that, lets have a look at the sizes of the vertices
in the graph to have something to compare to.
"""
@test [nout(start)] == nin(split) == [3nout(split)] == [sum(nin(joined))] == [nout(out)] == [6]
@test [nout(start), nout(joined)] == nin(out) == [6, 6]

# In many cases it is useful to hold on to the old graph before mutating
parentgraph = deepcopy(graph)

# It is not possible to change the size of `out` by exactly 2 due to `1.` and `2.` above.
# By default, NaiveNASlib warns when this happens and then tries to make the closest possible change.
# If we don't want the warning, we can tell NaiveNASlib to relax and make the closest possible change right away:
@test Δnout!(out => relaxed(2))

@test [nout(start)] == nin(split) == [3nout(split)] == [sum(nin(joined))] == [nout(out)] == [9]
@test [nout(start), nout(joined)] == nin(out) == [9, 9]

# As we can see above, the size change rippled through the graph due to the size relations described above.
# Pretty much every vertex was affected.
#
# Lets evaluate the graph just to verify that we don't get a dimension mismatch error.
@test graph((ones(6))) == [78, 78, 0, 114, 114, 0, 186, 186, 0]
# `TinyNNlib` uses [`parselect`](@ref) which will insert zeros when size increases by default.
# This helps the graph maintain the same function after mutation.
# In this case we changed the size of the output layer so we don't have
# the exact same function though, but hopefully it is clear why e.g. a
# linear layer after `out` would have made it produce the same output.
#
# Copy is still intact of course.
@test parentgraph((ones(6))) == [78, 78, 114, 114, 186, 186]

md"""
While we still have the complex model in scope, lets show a few more way to change the sizes. See the built in documentation for more information.

It is possible to supply a utility function for telling the utility of each neuron in a vertex.
NaiveNASlib will prioritize selecting the indices with higher utility.
"""
#
# Prefer high indices:
graphhigh = deepcopy(graph)
@test Δnout!(v -> 1:nout(v), graphhigh[end] => -3)
@test graphhigh((ones(6))) == [42, 0, 60, 0, 96, 0]

# Perfer low indices:
graphlow = deepcopy(graph)
@test Δnout!(v -> nout(v):-1:1, graphlow[end] => -3)
@test graphlow((ones(6))) == [78, 78, 114, 114, 186, 186]

# A simple approach when doing structured pruning is to prefer neurons with high magnitude.
# Here is how to set that as the default for LinearLayer.
# This is something one should probably implement in `TinyNNlib` instead.
using Statistics: mean
NaiveNASlib.defaultutility(l::LinearLayer) = mean(abs, l.W, dims=2)

graphhighmag = deepcopy(graph)
@test Δnout!(graphhighmag[end] => -3)
@test graphhighmag((ones(6))) == [78, 78, 114, 114, 186, 186]

# In many NAS applications one wants to apply random mutations to the graph.
# When doing so, one might end up in situations like this:
badgraphdecinc = deepcopy(graph)
v1, v2 = badgraphdecinc[[3, end]] # Imagine selecting these at random
@test Δnout!(v1 => relaxed(-2))
@test Δnout!(v2 => 6)
## Now we first deleted a bunch of weights, then we added new :(
@test badgraphdecinc((ones(6))) ==  [42, 0, 0, 60, 0, 0, 96, 0, 0]

# In such cases, it might be better to supply all wanted changes in one go and let
# NaiveNASlib try to come up with a decent compromise.
goodgraphdecinc = deepcopy(graph)
v1, v2 = goodgraphdecinc[[3, end]]
@test Δnout!(v1 => relaxed(-2), v2 => 3) # Mix relaxed and exact size changes freely
@test goodgraphdecinc((ones(6))) == [78, 78, 6, 0, 108, 114, 6, 6, 180, 180, 0, 0]

# It is also possible to change the input direction, but it requires specifying a size change for each input
# and is generally not recommended due to this.
graphΔnin = deepcopy(graph)
v1, v2 = graphΔnin[end-1:end]
## Use missing to signal "don't care"
@test Δnin!(v1 => (3, relaxed(2), missing), v2 => relaxed((1,2)))
@test nin(v1) == [6, 6, 6] # Sizes are tied to nout of split so they all have to be equal
@test nin(v2) == [18, 18] # Sizes are tied due to elementwise addition

# A common pruning strategy is to just remove the x% of params with lowest utility.
# This can be done by just not putting any size requirements and assign negative utility.
graphprune40 = deepcopy(graph)
Δsize!(graphprune40) do v
    utility = NaiveNASlib.defaultutility(v)
    ## We make some strong assumptions on weight distribution here for breviety :)
    return utility .- 0.4mean(utility)
end
@test nout.(vertices(graphprune40)) == [6, 6, 2, 2, 2, 2, 6, 6]
## Compare to original:
@test nout.(vertices(graph))        == [6, 9, 3, 3, 3, 3, 9, 9]

end #src

@testset "Weight modification example" begin #src

# ## Closer look at how weights are modified
# Here we take a closer look at how the weight matrices are changed. We use the following function to create
# `LinearLayer`s with easy to distinguish weights. It returns both a vertex and the `LinearLayer` so we
# can easily look at the weights:

function vertexandlayer(in, outsize)
    nparam = nout(in) * outsize
    l = LinearLayer(collect(reshape(1:nparam, :, nout(in))))
    return absorbvertex(l, in), l
end

# Make a simple model:
invertices = inputvertex.(["in1", "in2"], [3,4])
v1, l1 = vertexandlayer(invertices[1], 4)
v2, l2 = vertexandlayer(invertices[2], 3)
merged = conc(v1, v2, dims=1)
v3, l3 = vertexandlayer(invertices[1], nout(merged))
add = v3 + merged
v4, l4 = vertexandlayer(merged, 2)

# These weights might look a bit odd, but here we only care about making
# it easy to spot what has changed after size change below.
@test l1.W ==
[ 1 5  9
  2 6 10
  3 7 11
  4 8 12 ]

@test l2.W ==
[ 1 4 7 10
  2 5 8 11
  3 6 9 12 ]

@test l3.W ==
[ 1  8 15
  2  9 16
  3 10 17
  4 11 18
  5 12 19
  6 13 20
  7 14 21 ]

@test l4.W ==
[ 1 3 5 7  9 11 13 ;
  2 4 6 8 10 12 14 ]

# Now, lets decrease `v2` by 1 and force `merged` to retain its size
# which in turn forces `v1` to grow by 1.
# Assign utility 10 to neurons 1 and 3 of `v2` and 1 for all other neurons in the model.
@test Δnout!(v2 => -1, merged => 0) do v
    v == v2 ? [10, 1, 10] : 1
end

# `v1` got a new row of parameters at the end:
@test l1.W ==
[ 1 5  9
  2 6 10
  3 7 11
  4 8 12
  0 0  0 ]

# `v2` chose to drop its middle row as it was the output neuron with lowest utility:
@test l2.W ==
[ 1 4 7 10
  3 6 9 12 ]

# `v4` dropped the second to last column (which is aligned to the middle row of `v2`).
# and got new parameters in column 5 (which is aligned to the last row of `v1`):
@test l4.W ==
[  1 3 5 7 0  9 13
   2 4 6 8 0 10 14 ]

# Oh, and since `v4` is connected to `v3` though `add`, `v3` had the equivalent
# change to its rows:
@test l3.W ==
[ 1  8 15
  2  9 16
  3 10 17
  4 11 18
  0  0  0
  5 12 19
  7 14 21 ]
end #src

# ## Other modifications
# Lets just do a few quick examples of the other types of modifications.

# ### Add a vertex

# Using [`insert!`](@ref):

@testset "Add vertex example" begin #src
invertex = inputvertex("input", 3)
layer1 = linearvertex(invertex, 5)
graph = CompGraph(invertex, layer1)

## nvertices(g) is shortcut for length(vertices(g))
@test nvertices(graph) == 2
@test graph(ones(3)) == [3,3,3,3,3]

## Insert a layer between invertex and layer1
@test insert!(invertex, vertex -> linearvertex(vertex, nout(vertex))) # True if success

@test nvertices(graph) == 3
@test graph(ones(3)) == [9, 9, 9, 9, 9]
end #src

# ### Remove a vertex

# Using [`remove!`](@ref):
@testset "Remove vertex example" begin #src
invertex = inputvertex("input", 3)
layer1 = linearvertex(invertex, 5)
layer2 = linearvertex(layer1, 4)
graph = CompGraph(invertex, layer2)

@test nvertices(graph) == 3
@test graph(ones(3)) == [15, 15, 15, 15]

## Remove layer1 and change nin of layer2 from 5 to 3
## Would perhaps have been better to increase nout of invertex, but it is immutable
@test remove!(layer1) # True if success

@test nvertices(graph) == 2
@test graph(ones(3)) == [3, 3, 3, 3]
end #src

# ### Add an edge

# Using [`create_edge!`](@ref):
@testset "Add edge example" begin
invertices = inputvertex.(["input1", "input2"], [3, 2])
layer1 = linearvertex(invertices[1], 4)
layer2 = linearvertex(invertices[2], 4)
add = layer1 + layer2
out = linearvertex(add, 5)
graph = CompGraph(invertices, out)

@test nin(add) == [4, 4]
## Two inputs to this graph!
@test graph(ones(3), ones(2)) == [20, 20, 20, 20, 20]

## This graph is not interesting enough for there to be a good showcase for adding a new edge.
## Lets create a new layer which has a different output size just to see how things change
## The only vertex which support more than one input is add
layer3 = linearvertex(invertices[2], 6)
@test create_edge!(layer3, add) # True if success

## NaiveNASlib will try to increase the size in case of a mismatch by default
@test nin(add) == [6, 6, 6]
@test graph(ones(3), ones(2)) == [28, 28, 28, 28, 28]
end #src

# ### Remove an edge

# Using [`remove_edge!`](@ref):
@testset "Remove edge example" begin #src
invertex = inputvertex("input", 4)
layer1 = linearvertex(invertex, 3)
layer2 = linearvertex(invertex, 5)
merged = conc(layer1, layer2, layer1, dims=1)
out = linearvertex(merged, 3)
graph = CompGraph(invertex, out)

@test nin(merged) == [3, 5, 3]
@test graph(ones(4)) == [44, 44, 44]

@test remove_edge!(layer1, merged) # True if success

@test nin(merged) == [5, 3]
@test graph(ones(4)) == [32, 32, 32]
end #src
#

ext: jl


text:
using Statistics

export compute_hypercoef

"""
    compute_hypercoef( model )

compute hyperviscosity coefficient
"""
function compute_hypercoef(model)

    # compute gradient of w using 'gradient' matlab function

    # dxUx,dyUx = gradient(w(:,:,1)',grid.x_ref,grid.y_ref)
    # dxUy,dyUy = gradient(w(:,:,2)',grid.x_ref,grid.y_ref)

    nx, ny = model.grid.nx, model.grid.ny

    update_velocities!(model)

    dxUx = irfft(1im .* model.grid.kx .* model.û_x, nx)
    dyUx = irfft(1im .* model.grid.ky .* model.û_x, nx)
    dxUy = irfft(1im .* model.grid.kx .* model.û_y, nx)
    dyUy = irfft(1im .* model.grid.ky .* model.û_y, nx)

    s = (dxUx .+ dyUy) ./ 2
    d = (dyUx .+ dxUy) ./ 2

    @show lambda_rms = sqrt(mean(s .^ 2 .+ d .^ 2))

    dx, dy = model.grid.dx, model.grid.dy

    # Hyperviscosity coefficient
    coef = 40 * lambda_rms * ((dx+dy)/2π)^model.hv_order

    return coef

end


text:
using IcanProj
using JLD2

###################################################################################

# Settings

opt = get_options()

###################################################################################

# Load metadata about watersheds

@load Pkg.dir("IcanProj", "data", "wsh_info.jld2") wsh_info

# Read soil parameter file

soil_param = read_soil_params(opt["base_folder"])

# Read vegetation parameter file

veg_param = read_veg_param(opt["base_folder"])

###################################################################################

# Clear target folder

rm(opt["target_folder"], force = true, recursive = true)

###################################################################################

# Loop over selected watersheds

for wsh_single in wsh_info

    # Loop over grid resolutions

    for res in [:ind_1km, :ind_5km, :ind_10km, :ind_25km, :ind_50km]

        try

            info("Processing $(wsh_single.name) for resolution $(string(res)[5:end])")

            # Folder names

            path_sim = joinpath(opt["target_folder"], string(wsh_single.regine_main), string(res)[5:end])

            path_params = joinpath(path_sim, "params")

            path_results = joinpath(path_sim, "results")

            path_forcings = joinpath(path_sim, "forcing")

            # Create folders

            mkpath(path_params)

            mkpath(path_results)

            mkpath(path_forcings)

            # Write global parameter file

            if res == :ind_1km

                write_global_param(path_sim,
                                   opt["base_folder"],
                                   opt["startyear"],
                                   opt["endyear"],
                                   opt["timestep"],
                                   opt["output_force"],
                                   opt["full_energy"],
                                   opt["output_binary"])

            else

                write_global_param(path_sim,
                                   path_sim,
                                   opt["startyear"],
                                   opt["endyear"],
                                   opt["timestep"],
                                   opt["output_force"],
                                   opt["full_energy"],
                                   opt["output_binary"])

            end

            # Write vegetation library file

            file_veg_src = joinpath(opt["base_folder"], "params", "veglib_param")
            file_veg_dst = joinpath(path_params, "veglib_param")

            cp(file_veg_src, file_veg_dst, remove_destination=true)

            # Aggregate and write soil parameter file

            soil_param_ave = average_soilparams(soil_param, wsh_single, res)

            write_soil_params(path_sim, soil_param_ave)

            # Write vegetation parameter file

            veg_param_ave = average_vegparams(veg_param, wsh_single, res)

            write_veg_params(path_sim, veg_param_ave)

            # Write model forcing data

            if res != :ind_1km

                average_forcings(opt["base_folder"], path_sim, soil_param, wsh_single, res)

            end

        catch

            info("Failed processing $(wsh_single.name) for resolution $(string(res)[5:end])")

        end

    end

end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

#   hadamard_sum!(du, A, volume_flux, orientation, u, equations,
#                              sparsity_pattern, skip_index=(i,j)->false)
#
# Computes the flux difference ∑_j A[i, j] * f(u_i, u_j) and accumulates the result into `du`.
# - `du`, `u` are vectors
# - `A` is the skew-symmetric flux differencing matrix.
# - `sparsity_pattern`: either `nothing` or an AbstractSparseMatrix which specifies the sparsity
#   pattern of `A`
@inline function hadamard_sum!(du, A, volume_flux, orientation, u, equations,
                               sparsity_pattern::Nothing, skip_index=(i,j)->false)

  rows, cols = axes(A)
  for i in rows
    u_i = u[i]
    du_i = du[i]
    for j in cols
      # This routine computes only the upper-triangular part of the hadamard sum (A .* F).
      # We avoid computing the lower-triangular part, and instead accumulate those contributions
      # while computing the upper-triangular part (using the fact that A is skew-symmetric and F
      # is symmetric).
      if j > i && !skip_index(i, j)
          AF_ij = A[i,j] * volume_flux(u_i, u[j], orientation, equations)
          du_i = du_i + AF_ij
          du[j] = du[j] - AF_ij
      end
    end
    du[i] = du_i
  end
end

@inline function hadamard_sum!(du, A, volume_flux, orientation, u,
                               equations, sparsity_pattern::AbstractSparseMatrix{Bool})
  n = size(sparsity_pattern, 2)
  rows = rowvals(sparsity_pattern)
  for i in 1:n
    u_i = u[i]
    du_i = du[i]
    for id in nzrange(sparsity_pattern, i)
      j = rows[id]
      # This routine computes only the upper-triangular part of the hadamard sum (A .* F).
      # We avoid computing the lower-triangular part, and instead accumulate those contributions
      # while computing the upper-triangular part (using the fact that A is skew-symmetric and F
      # is symmetric).
      if j > i
        AF_ij = A[i,j] * volume_flux(u_i, u[j], orientation, equations)
        du_i = du_i + AF_ij
        du[j] = du[j] - AF_ij
      end
    end
    du[i] = du_i
  end
end

# For DGMulti implementations, we construct "physical" differentiation operators by taking linear
# combinations of reference differentiation operators scaled by geometric change of variables terms.
# We use a lazy evaluation of physical differentiation operators, so that we can compute linear
# combinations of differentiation operators on-the-fly in an allocation-free manner.
function build_lazy_physical_derivative(element, orientation,
                                        mesh::VertexMappedMesh{2}, dg, cache)
  @unpack Qrst_skew = cache
  @unpack rxJ, sxJ, ryJ, syJ = mesh.md
  if orientation == 1
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (rxJ[1,element], sxJ[1,element]))
  else # if orientation == 2
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (ryJ[1,element], syJ[1,element]))
  end
end

function build_lazy_physical_derivative(element, orientation,
                                        mesh::VertexMappedMesh{3}, dg, cache)
  @unpack Qrst_skew = cache
  @unpack rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ = mesh.md
  if orientation == 1
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (rxJ[1,element], sxJ[1,element], txJ[1,element]))
  elseif orientation == 2
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (ryJ[1,element], syJ[1,element], tyJ[1,element]))
  else # if orientation == 3
    return LazyMatrixLinearCombo(Qrst_skew, 2 .* (rzJ[1,element], szJ[1,element], tzJ[1,element]))
  end
end

function compute_flux_differencing_SBP_matrices(dg::DGMulti)
  rd = dg.basis
  Qrst_hybridized, VhP, Ph = StartUpDG.hybridized_SBP_operators(rd)
  Qrst_skew = map(A -> 0.5 * (A - A'), Qrst_hybridized)
  return Qrst_skew, VhP, Ph
end

function compute_flux_differencing_SBP_matrices(dg::DGMultiFluxDiff{<:SBP})
  rd = dg.basis
  @unpack M, Drst, Pq = rd
  Qrst = map(D -> M * D, Drst)
  Qrst_skew = map(A -> 0.5 * (A - A'), Qrst)
  return Qrst_skew
end

# precompute sparsity pattern for optimized flux differencing routines for tensor product elements
function compute_sparsity_pattern(flux_diff_matrices, dg::DG,
                                  tol = 100 * eps(real(dg))) where {DG <: DGMultiFluxDiff{<:SBP, <:Union{Quad, Hex}}}
  sparsity_pattern = sum(map(A->abs.(A), droptol!.(sparse.(flux_diff_matrices), tol))) .!= 0
  return sparsity_pattern
end

compute_sparsity_pattern(flux_diff_matrices, dg::DGMulti) = nothing

# For flux differencing SBP-type approximations, store solutions in Matrix{SVector{nvars}}.
# This results in a slight speedup for `calc_volume_integral!`.
function allocate_nested_array(uEltype, nvars, array_dimensions, dg::DGMultiFluxDiff{<:SBP})
  return zeros(SVector{nvars, uEltype}, array_dimensions...)
end

# Initialization specialized to solution storage for DGMultiFluxDiff with SBP approximation type.
function compute_coefficients!(u, initial_condition, t,
                               mesh::AbstractMeshData{NDIMS}, equations, dg::DGMultiFluxDiff{<:SBP}, cache) where {NDIMS}
  md = mesh.md
  rd = dg.basis

  # evaluate the initial condition at nodal point
  @threaded for i in each_quad_node_global(mesh, dg, cache)
    u[i] = initial_condition(getindex.(md.xyzq, i), t, equations)
  end
end

function create_cache(mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:SBP}, RealT, uEltype)

  rd = dg.basis
  md = mesh.md

  # for use with flux differencing schemes
  Qrst_skew = compute_flux_differencing_SBP_matrices(dg)
  sparsity_pattern = compute_sparsity_pattern(Qrst_skew, dg)

  nvars = nvariables(equations)

  # Todo: simplices. Factor common storage into a struct (MeshDataCache?) for reuse across solvers?
  # storage for volume quadrature values, face quadrature values, flux values
  u_values = allocate_nested_array(uEltype, nvars, size(md.xq), dg)
  u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
  flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
  lift_scalings = rd.wf ./ rd.wq[rd.Fmask] # lift scalings for diag-norm SBP operators

  local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) for _ in 1:Threads.nthreads()]

  # Use an array of SVectors (chunks of `nvars` are contiguous in memory) to speed up flux differencing
  fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, rd.Nq) for _ in 1:Threads.nthreads()]

  return (; md, Qrst_skew, sparsity_pattern, invJ = inv.(md.J), lift_scalings, inv_wq = inv.(rd.wq),
            u_values, u_face_values, flux_face_values,
            local_values_threaded, fluxdiff_local_threaded)
end

function create_cache(mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:Polynomial}, RealT, uEltype)

  rd = dg.basis
  @unpack md = mesh

  Qrst_skew, VhP, Ph = compute_flux_differencing_SBP_matrices(dg)
  sparsity_pattern = compute_sparsity_pattern(Qrst_skew, dg)

  nvars = nvariables(equations)

  # temp storage for entropy variables at volume quad points
  entropy_var_values = allocate_nested_array(uEltype, nvars, (rd.Nq, md.num_elements), dg)

  # storage for all quadrature points (concatenated volume / face quadrature points)
  num_quad_points_total = rd.Nq + rd.Nfq
  entropy_projected_u_values = allocate_nested_array(uEltype, nvars, (num_quad_points_total, md.num_elements), dg)
  projected_entropy_var_values = allocate_nested_array(uEltype, nvars, (num_quad_points_total, md.num_elements), dg)

  # initialize temporary storage as views into entropy_projected_u_values
  u_values = view(entropy_projected_u_values, 1:rd.Nq, :)
  u_face_values = view(entropy_projected_u_values, rd.Nq+1:num_quad_points_total, :)
  flux_face_values = similar(u_face_values)

  # local storage for interface fluxes, rhs, and source
  local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) for _ in 1:Threads.nthreads()]

  # Use an array of SVectors (chunks of `nvars` are contiguous in memory) to speed up flux differencing
  # The result is then transferred to rhs_local_threaded::StructArray{<:SVector} before
  # projecting it and storing it into `du`.
  fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, num_quad_points_total) for _ in 1:Threads.nthreads()]
  rhs_local_threaded = [allocate_nested_array(uEltype, nvars, (num_quad_points_total,), dg)  for _ in 1:Threads.nthreads()]

  return (; md, Qrst_skew, sparsity_pattern, VhP, Ph, invJ = inv.(md.J),
            entropy_var_values, projected_entropy_var_values, entropy_projected_u_values,
            u_values, u_face_values,  flux_face_values,
            local_values_threaded, fluxdiff_local_threaded, rhs_local_threaded)
end

function entropy_projection!(cache, u, mesh::VertexMappedMesh, equations, dg::DGMulti)

  rd = dg.basis
  @unpack Vq = rd
  @unpack VhP, entropy_var_values, u_values, entropy_var_values = cache
  @unpack projected_entropy_var_values, entropy_projected_u_values = cache

  # TODO: simplices. Address hard-coding of `entropy2cons!`
  apply_to_each_field(mul_by!(Vq), u_values, u)
  entropy_var_values .= cons2entropy.(u_values, equations)

  # "VhP" fuses the projection "P" with interpolation to volume and face quadrature "Vh"
  apply_to_each_field(mul_by!(VhP), projected_entropy_var_values, entropy_var_values)
  entropy_projected_u_values .= entropy2cons.(projected_entropy_var_values, equations)
end

function calc_volume_integral!(du, u, volume_integral,
                               mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:SBP}, cache)

  rd = dg.basis
  @unpack fluxdiff_local_threaded, sparsity_pattern, inv_wq = cache
  @unpack volume_flux = volume_integral

  # Todo: simplices. Dispatch on curved/non-curved mesh types, this code only works for affine meshes (accessing rxJ[1,e],...)
  @threaded for e in eachelement(mesh, dg, cache)
    fluxdiff_local = fluxdiff_local_threaded[Threads.threadid()]
    fill!(fluxdiff_local, zero(eltype(fluxdiff_local)))
    u_local = view(u, :, e)
    for i in eachdim(mesh)
      Qi_skew = build_lazy_physical_derivative(e, i, mesh, dg, cache)

      # use `sparsity_pattern` to dispatch hadamard_sum!. If using Tri or Tet elements, the
      # Qi_skew matrices are dense, and sparsity_pattern = nothing. If using Quad or Hex
      # elements with an SBP approximationType, then sparsity_pattern::AbstractSparseMatrix{Bool}.
      hadamard_sum!(fluxdiff_local, Qi_skew, volume_flux, i,
                    u_local, equations, sparsity_pattern)
    end

    for i in each_quad_node(mesh, dg, cache)
      du[i, e] = du[i, e] + fluxdiff_local[i] * inv_wq[i]
    end
  end
end

function calc_volume_integral!(du, u, volume_integral,
                               mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:Polynomial}, cache)

  rd = dg.basis
  @unpack entropy_projected_u_values, Ph, sparsity_pattern = cache
  @unpack fluxdiff_local_threaded, rhs_local_threaded = cache
  @unpack volume_flux = volume_integral

  # skips subblock of Qi_skew which we know is zero by construction
  skip_index(i,j) = i > rd.Nq && j > rd.Nq

  # Todo: simplices. Dispatch on curved/non-curved mesh types, this code only works for affine meshes (accessing rxJ[1,e],...)
  @threaded for e in eachelement(mesh, dg, cache)
    fluxdiff_local = fluxdiff_local_threaded[Threads.threadid()]
    fill!(fluxdiff_local, zero(eltype(fluxdiff_local)))
    u_local = view(entropy_projected_u_values, :, e)
    for i in eachdim(mesh)
      Qi_skew = build_lazy_physical_derivative(e, i, mesh, dg, cache)
      hadamard_sum!(fluxdiff_local, Qi_skew, volume_flux, i,
                    u_local, equations, sparsity_pattern, skip_index)
    end

    # convert fluxdiff_local::Vector{<:SVector} to StructArray{<:SVector}
    rhs_local = rhs_local_threaded[Threads.threadid()]
    for i in Base.OneTo(length(fluxdiff_local))
      rhs_local[i] = fluxdiff_local[i]
    end
    apply_to_each_field(mul_by_accum!(Ph), view(du, :, e), rhs_local)
  end
end

# Specializes on Polynomial (e.g., modal) DG methods with a flux differencing volume kernel, e.g.,
# an entropy conservative/stable discretization. For modal DG schemes, an extra `entropy_projection`
# is required (see https://doi.org/10.1016/j.jcp.2018.02.033, Section 4.3).
function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions::BC,
              source_terms::Source, dg::DGMultiFluxDiff{Polynomial}, cache) where {Source, BC}

  @trixi_timeit timer() "Reset du/dt" fill!(du, zero(eltype(du)))

  # this function evaluates the solution at volume and face quadrature points (which was previously
  # done in `prolong2interfaces` and `calc_volume_integral`)
  @trixi_timeit timer() "entropy_projection!" entropy_projection!(cache, u, mesh, equations, dg)

  @trixi_timeit timer() "calc_volume_integral!" calc_volume_integral!(du, u, dg.volume_integral,
                                                                      mesh, equations, dg, cache)

  # the following functions are the same as in VolumeIntegralWeakForm, and can be reused from dg.jl
  @trixi_timeit timer() "calc_interface_flux!" calc_interface_flux!(cache, dg.surface_integral,
                                                                    mesh, equations, dg)

  @trixi_timeit timer() "calc_boundary_flux!" calc_boundary_flux!(cache, t, boundary_conditions,
                                                                  mesh, equations, dg)

  @trixi_timeit timer() "calc_surface_integral!" calc_surface_integral!(du, u, dg.surface_integral,
                                                                        mesh, equations, dg, cache)

  @trixi_timeit timer() "invert_jacobian" invert_jacobian!(du, mesh, equations, dg, cache)

  @trixi_timeit timer() "calc_sources!" calc_sources!(du, u, t, source_terms,
                                                      mesh, equations, dg, cache)

  return nothing
end

end # @muladd


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

mutable struct IoK8sApiCoreV1ConfigMapEnvSource <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    optional::Any # spec type: Union{ Nothing, Bool } # spec name: optional

    function IoK8sApiCoreV1ConfigMapEnvSource(;name=nothing, optional=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ConfigMapEnvSource, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1ConfigMapEnvSource, Symbol("optional"), optional)
        setfield!(o, Symbol("optional"), optional)
        o
    end
end # type IoK8sApiCoreV1ConfigMapEnvSource

const _property_map_IoK8sApiCoreV1ConfigMapEnvSource = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("optional")=>Symbol("optional"))
const _property_types_IoK8sApiCoreV1ConfigMapEnvSource = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("optional")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1ConfigMapEnvSource }) = collect(keys(_property_map_IoK8sApiCoreV1ConfigMapEnvSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1ConfigMapEnvSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1ConfigMapEnvSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ConfigMapEnvSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ConfigMapEnvSource[property_name]

function check_required(o::IoK8sApiCoreV1ConfigMapEnvSource)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ConfigMapEnvSource }, name::Symbol, val)
end


text:
function arterialpressure!(system::CVSystem,n::Int64)
    for i = 1:length(system.branches.ID)
        system.branches.P[i][:,n+2] .= system.branches.beta[i][end].*
            (system.branches.A[i][:,n+2].^0.5.-
            sqrt.(system.branches.A0[i][end]))./mmHgToPa;
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
function jump_max_area_perimeter_Dn(n)
    if mod(log2(n),1)==0 && n>=16
        d = atan(tan(2*pi/n)*tan(pi/n))-asin(sin(2*pi/n)^2/sqrt(2+2*cos(2*pi/n)*cos(4*pi/n)));
        Ll = 2*n*sin(pi/(2*n))*cos(d/2);
        Lu = 2*n*sin(pi/(2*n));
        c = zeros(Int(3*n/8));
        for i = 1:Int(3*n/8)
            c[i] = 1 + round(mod(i,3)/3);
        end
        model = Model(Ipopt.Optimizer)
        #model = Model(() -> AmplNLWriter.Optimizer(Couenne_jll.amplexe))
        # Variables
        @variable(model, a[i = 1:Int(3*n/8)] >= 0, start = pi/n+(-1)^(i-1)*d)
        # Perimeter
        @NLobjective(model, Max, sum(4*c[i]*sin(a[i]/2) for i in 1:Int(3*n/8)))
        # Optimal diameter graph
        @constraint(model, sum(c[i]*a[i] for i in 1:Int(3*n/8)) == pi/2)
        @NLconstraint(model, sum((-1)^(i-1)*sin(sum(c[j]*a[j] for j in 1:i)) for i in 1:Int(3*n/8-1)) == 1/2)
        @constraint(model, [i = 1:Int(3*n/8)], c[i]*a[i] <= pi/3)
        # Bounds
        @NLconstraint(model, Ll <= sum(4*c[i]*sin(a[i]/2) for i in 1:Int(3*n/8)) <= Lu)
        # Solving
        optimize!(model)
        L = objective_value(model)
        a = value.(a)
        ct = solve_time(model)
        return(L,a,ct)
    end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
if VERSION < v"0.7.0" && Pkg.installed("BinDeps") == nothing
    Pkg.add("BinDeps")
end

using BinDeps
using Compat.Libdl

@BinDeps.setup

version = "9.1.7"
bilevelBenchmark = library_dependency("bilevelBenchmarkJulia",
                                      aliases=["blb18_op_v$version", "blb18_op.exe"],
                                      os = BinDeps.OSNAME)

# build from source
provides(Sources,
        URI("https://github.com/jmejia8/bilevel-benchmark/archive/v$(version).zip"),
        unpacked_dir="bilevel-benchmark-$version",
        bilevelBenchmark, so = :Unix)

prefix = joinpath(BinDeps.depsdir(bilevelBenchmark), "usr")
srcdir = joinpath(BinDeps.srcdir(bilevelBenchmark), "bilevel-benchmark-$version")

provides(SimpleBuild,
    (@build_steps begin
        GetSources(bilevelBenchmark)
        CreateDirectory(joinpath(prefix, "lib"))
        @build_steps begin
            ChangeDirectory(srcdir)
            MAKE_CMD
            `mv blb18_op.$(Libdl.dlext) "$prefix/lib/blb18_op_v$(version).$(Libdl.dlext)"`
        end
    end), [bilevelBenchmark], os = :Unix)

provides(Binaries,
    URI("https://github.com/jmejia8/bilevel-benchmark/files/5967963/blb_op.zip"),
    bilevelBenchmark,
    os = :Windows)

@BinDeps.install Dict(:bilevelBenchmarkJulia => :bilevelBenchmark)


text:
#
## Purity heuristic
#
#  Next to the PPT condition we use a purity enforcing heuristic to find choi states corresponding
#  to RO operation in which the size of the environment is known.
#
# Input:
#    rho   - the input state rhoAB to be distilled, assuming dim(A)=dim(B)
#    dout  - desired dimension of the maximally entangled output state (on one node).
#    nE_A  - size of the environment to be used at node A
#    nE_B  - size of the environment to be used at node B
#    maxIter - maximum number of iterations
#

using Convex
using SCS

export purity

function purity(rho::AbstractArray, dout::Int, nE_A::Int, nE_B::Int, maxIter::Int)

	# Define a tolerance parameter: eigenvalues below that number are considered 0
	tol = (10//1)^(-5);

	# In analogy to the notation used in the paper, set the output dimensions of Ahat and Bhat
        # to contain the (ideally) maximally entangled output state.
  	nAhat = dout;
  	nBhat = nAhat;

	# Dimension of the input state
	(d1,d2) = size(rho);
	din = convert(Int,sqrt(d1));
  	nAprime = din;
  	nBprime = nAprime;

	# Dimension of the flag registers, only 2 possible values: success or failure
  	nF_A = 2;
  	nF_B = 2;

	# Compute the total dimension of the A and B maps
  	nA = convert(Int,nAhat * nF_A * nAprime * nE_A);
  	nB = convert(Int,nBhat * nF_B * nBprime * nE_B);
  	nC = nA * nB;

	# vector of system dimensions: specifies how we order them
	sysDims = [nAhat, nBhat, nF_A, nF_B, nE_A, nE_B, nAprime, nBprime];

	# At the first level of the iteration, no additional constraints will be placed
  	t_A = 1;
  	t_B = 1;

  	# instantiate direction matrices
  	W_Abar = zeros(nA, nA);
  	W_Bbar = zeros(nB, nB);

  	# do iterations until purity is reached or 100 iterations were done
  	for i = 1:maxIter
    		println("Starting iteration ", i)

		# Solve the optimization problem
    		@time (C_AbarBbar, F, p_succ) = findChoi(rho, W_Abar, W_Bbar, t_A, t_B, nE_A, nE_B, dout)

    		# sometimes Convex.jl fails.
    		# try again in that case
    		j = 1;
    		while C_AbarBbar == zeros((nC, nC))
      			println(" trying again ", j, "...")
      			j += 1
      			@time (C_AbarBbar, F, p_succ) = findChoi(rho, W_Abar, W_Bbar, t_A, t_B, nE_A, nE_B, dout);
    		end

		# Compute the marginal states. Ideally these are pure.
		C_Abar = partialtrace(CABarBbar,  [2, 4, 6, 8], sysDims);
		C_Bbar = partialtrace(CABarBbar,  [1, 3, 5, 7], sysDims);

		# Verify whether the A and B part are pure
    		eA = eigvals(C_Abar);
    		eA = -sort(-real(eA));
    		println(" eA = ", round(eA[1:10], 2));

		eB = eigvals(C_Bbar);
		eB = -sort(-real(eB));
		println(" eB = ", round(eB[1:10],2));

		# Print the fidelity and success probability attained
    		println(" F = ", F)
    		println(" p_succ = ", p_succ)

    		# check with some tolerance if the state is pure
    		# if so, return the found Choi state
    		if (eA <= tol) && (eB <= tol)
      			return (C_Abar,C_Bbar);
    		end

    		@time (W_Abar, t_A) = findW(C_Abar);
    		@time (W_Bbar, t_B) = findW(C_Bbar);

  	end

  	return(C_Abar, C_Bbar);
end

#
# Solve the original SDP, including the additional constraint on the matrices enforcing the purity heuristic
#
# Input:
#   rho		- input state to be distilled
#   W_Abar	- matrix enforcing the purity heuristic on system A
#   W_Bbar      - matrix enforcing the purity heuristic on system B
#   t_A		- constraint t_A enforcing purity on system A
#   t_B		- constraint t_B enforcing purity on system B
#   nE_A 	- size of the environment on A
#   nE_B 	- size of the environment on B
#   dout	- dimension of the output system (on one side)
#
# Output:
#
#   C		- overall choi states
#   F		- fidelity achieved
#   psucc	- success probability achieved
#

function findChoi(rho::AbstractArray, W_Abar::AbstractArray, W_Bbar::AbstractArray, t_A::Number, t_B::Number, nE_A::Number, nE_B::Number, dout::Number)

	# verify the input
	@assert isQuantumState(rho) "expected rho to be a quantum state"
  	@assert 0 <= t_A <= 1
  	@assert 0 <= t_B <= 1

	# In analogy to the notation used in the paper, set the output dimensions of Ahat and Bhat
        # to contain the (ideally) maximally entangled output state.
  	nAhat = dout;
  	nBhat = nAhat;

	# Dimension of the input state
	(d1,d2) = size(rho);
	din = convert(Int,sqrt(d1));
  	nAprime = din;
  	nBprime = nAprime;

	# Dimension of the flag registers, only 2 possible values: success or failure
  	nF_A = 2;
  	nF_B = 2;

	# Compute the total dimension of the A and B maps
  	nA = convert(Int,nAhat * nF_A * nAprime * nE_A);
  	nB = convert(Int,nBhat * nF_B * nBprime * nE_B);
  	nC = convert(Int,nA * nB);

	# vector of system dimensions: specifies how we order them
	sysDims = [nAhat, nBhat, nF_A, nF_B, nE_A, nE_B, nAprime, nBprime];

	# At the first level of the iteration, no additional constraints will be placed
  	t_A = 1;
  	t_B = 1;

	# Compute projector onto the 11 flag indicating success
  	e0 = eVec(2, 1);
  	e1 = eVec(2, 2);

	e11 = kron(e1,e1);
  	e11 = e11 * e11';

	# Compute a projector onto the Epr pair
	epr = maxEnt(nAhat);

	# Set up the optimization problem. This SDP optimizes (For now XXX) the product of the
	# fidelity and the success probability

	# The overall choi state
  	C_AbarBbar = Semidefinite(nC);

	# The partial trace components on A and B
  	C_Abar = partialtrace(C_AbarBbar, [2, 4, 6, 8], sysDims);
  	C_Bbar = partialtrace(C_AbarBbar, [1, 3, 5, 7], sysDims);

  	C_AprimeBprime = partialtrace(C_AbarBbar, [1, 2, 3, 4, 5, 6] , sysDims);

	# Define the objective function: maximize F * psucc
  	problem = maximize(trace(kron(kron(epr,eye(nE_A * nE_B)),kron(e11,rho')) * C_AbarBbar));

  	# constraints from purity
  	problem.constraints += [trace(C_Abar * W_Abar) <= t_A];
  	problem.constraints += [trace(C_Bbar * W_Bbar) <= t_B];

  	# constraints for Choi state
  	problem.constraints += [C_AprimeBprime == eye(din^2)/din^2];
  	problem.constraints += [trace(C_AbarBbar) == 1];

  	# PPT constraint on B
	print("Before PPT\n");
  	problem.constraints += partialtranspose(C_AbarBbar, [2, 4, 6, 8], sysDims) in :SDP;

	print("Running SDP\n");
  	solve!(problem, SCSSolver(verbose = true));

  	p_succ = nAprime * nBprime  * trace(kron(kron(eye(nAhat * nBhat * nE_A * nE_B),e11),rho') * CAbarBbar.value);
  	F = problem.optval * nAprime * nBprime/ p_succ;

  	return (C_AbarBbar.value, F, p_succ);
end

#
## findW
#
# Solves the iterative SDP of the purity heuristic, trying to force the relevant part of the input
# to be pure by searching for an updated W. The new t will be the new SDP value
#

function findW(C_opt::AbstractArray, verbose::Bool = false)

	# Make sure the input is actually a quantum state
  	@assert isQuantumState(C_opt) "expected C_opt to be a quantum state"

	# Recover the size of the input
  	nC = size(C_opt)[1];

	# Search for the new W
  	W_Xbar = Semidefinite(nC);

	# Define the objective
  	problem = minimize(trace(C_opt' * W_Xbar));

	# Normalize W
  	problem.constraints += ([W_Xbar <= eye(nC)]);

	# If W was pure itself, then it'll have trace 1
  	problem.constraints += ([trace(W_Xbar) == nC - 1]);

  	solve!(problem, SCSSolver(verbose = verbose));
  	return (W_Xbar.value, problem.optval)
end


text:
# Copyright (c) 2015: AmplNLWriter.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

using Test

@testset "MOI_wrapper" begin
    include("MOI_wrapper.jl")
end


text:
"""
    IWidthable(nexid)

An interface for widgets whom width can modified.
"""
struct IWidthable <: AbstractINextion
    _nid::NexID

    function IWidthable(nid::NexID)
        new(nid)
    end
end

"""
    obj.width = new_val

Set Nextion object `obj` width to `new_val`.
"""
function setproperty!(obj::IWidthable, property::Symbol, new_val)
    if property == :width
        setnexproperty!(NexID(obj), :wid, new_val)
    else
        setfield!(obj, property, new_val)
    end
end


text:
using Dash

app = dash(meta_tags = [Dict(["charset"=>"ISO-8859-1", "keywords"=>"dash,pleasant,productive", "http-equiv"=>"content-type", "content"=>"text/html"])])

app.layout = html_div(children = "Hello world!", id = "hello-div")

run_server(app)


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#
# ------------------------------------------
# mostly static tests of multithreading
#
sleeptime = 0.5
a = [0.0]
b = [0.0]
c = [0.0]
d = [0.0]
incr!(x) = x[1] += 1

println("... testing multithreading  stage 1, (sleeptime=$sleeptime) ...")
println("number of available threads: ", nthreads())

clk = PClock()
Δt = clk.Δt
sleep(sleeptime)

print(clk)
@test clk.id == 1
m = match(r"Clock 1 \(\+(\d+)\)", repr(clk))
@test parse(Int, m.captures[1]) == nthreads()-1
@test length(clk.ac) == nthreads()-1
@test clk.ac[1].thread == 2

# parallel clock identification
id = nthreads()+1
@test_warn "parallel clock $id not available!" pclock(clk, id)
c2 = pclock(clk, 2)
@test c2.id == 2
@test c2.clock.id == 2
@test c2.clock.ac[].master[].id == 1
m = match(r"Active clock (\d+)\:", repr(c2))
@test parse(Int, m.captures[1]) == 2
@test pclock(c2.clock, 2) isa DiscreteEvents.ActiveClock
@test clk == pclock(c2.clock, 1)

# test fork!
@test_warn "only the master clock" fork!(c2.clock)

# test _cid with parallel clocks
@test DiscreteEvents._cid(clk, 2, false) == 2
@test DiscreteEvents._cid(clk, 2, true) in 1:nthreads()
@test DiscreteEvents._cid(clk, nthreads()+1, false) == 1
@test DiscreteEvents._cid(clk, nthreads()+1, true) in 1:nthreads()
@test DiscreteEvents._cid(c2, 2, false) == 2
@test DiscreteEvents._cid(c2, 2, true) in 1:nthreads()
@test DiscreteEvents._cid(c2, 1, false) == 1
@test DiscreteEvents._cid(c2, 1, true) == 1
@test DiscreteEvents._cid(c2, nthreads()+1, false) == 2
@test DiscreteEvents._cid(c2, nthreads()+1, true) in 1:nthreads()
@test DiscreteEvents._cid(c2.clock, 2, false) == 2
@test DiscreteEvents._cid(c2.clock, 2, true) in 1:nthreads()
@test DiscreteEvents._cid(c2.clock, 1, false) == 1
@test DiscreteEvents._cid(c2.clock, 1, true) == 1
@test DiscreteEvents._cid(c2.clock, nthreads()+1, false) == 2
@test DiscreteEvents._cid(c2.clock, nthreads()+1, true) in 1:nthreads()

# exception handling and diagnosis
if DiscreteEvents._handle_exceptions[1]
    println("... remote error handling ...")
    put!(clk.ac[1].forth, DiscreteEvents.Clear())
    sleep(sleeptime)
    err = diagnose(clk, 2)
    @test err[1] isa ErrorException
    # @test occursin(r"^step\!\(\:\:DiscreteEvents.ActiveClock.+ at threads\.jl", string(err[2][2]))
end
old = DiscreteEvents._handle_exceptions[1]
DiscreteEvents._handle_exceptions[1] = false
put!(clk.ac[1].forth, DiscreteEvents.Clear())
sleep(sleeptime)
@test diagnose(clk, 2).state === :failed
DiscreteEvents._handle_exceptions[1] = old
@test_warn "parallel clock $id not available!" diagnose(clk, id)

# now testing parallel events and runs
#
# start a parallel clock again
clk = PClock()
Δt = clk.Δt
sleep(sleeptime)
c2 = pclock(clk, 2)

function pinc!(c::Clock, x)
    delay!(c, 1)
    incr!(x)
end

println("... register events to parallel clock 2 over master ...")
event!(clk, fun(incr!, a), 1.0, cid=2)
sleep(sleeptime)
@test length(c2.clock.sc.events) == 1
event!(clk, fun(incr!, b), fun(≥, fun(tau, c2), 5), cid=2)
sleep(sleeptime)
@test length(c2.clock.sc.cevents) == 1
periodic!(clk, fun(incr!, c), cid=2)
sleep(sleeptime)
@test length(c2.clock.sc.samples) == 1
process!(clk, Prc(1, pinc!, d), cid=2)
sleep(sleeptime)
@test length(c2.clock.processes) == 1

println("... 1st run parallel clock 2 (thread 2) ...")
# put!(clk.ac[1].forth, DiscreteEvents.Run(10.0))
# @test take!(clk.ac[1].back).t > 0
# sleep(sleeptime)
tns = [0]
t1 = time_ns()
iter = [0]
while (c2.clock.time+Δt) ≤ 10
    put!(clk.ac[1].forth, DiscreteEvents.Run(Δt, false))
    tns[1] += take!(clk.ac[1].back).t
    iter[1] += 1
end
t1 = time_ns() - t1
println("$iter ticks took $(Int(t1)*1e-9) s, clock time $(Int(tns[1])*1e-9) s")
@test c2.clock.time ≈ 10
@test c2.clock.scount == 1000
@test a[1] == 1
@test b[1] == 1
@test c[1] == 1000
@test d[1] >= 10

println("... parallel clock 2 reset ...")
put!(clk.ac[1].forth, DiscreteEvents.Reset(true))
@test take!(clk.ac[1].back).x == 1
@test c2.clock.time == 0
@test c2.clock.scount == 0
empty!(c2.clock.processes)

println("... register events to active clock 2 ...")
event!(c2, fun(incr!, a), 1.0)
sleep(sleeptime)
@test length(c2.clock.sc.events) == 1
event!(c2, fun(incr!, b), fun(≥, fun(tau, c2), 5))
sleep(sleeptime)
@test length(c2.clock.sc.cevents) == 1
periodic!(c2, fun(incr!, c))
sleep(sleeptime)
@test length(c2.clock.sc.samples) == 1

println("... 2nd run parallel clock 2 (thread 2) ...")
# put!(clk.ac[1].forth, DiscreteEvents.Run(10.0))
# @test take!(clk.ac[1].back).t > 0
# sleep(sleeptime)
tns = [0]
t1 = time_ns()
iter = [0]
while (c2.clock.time+Δt) ≤ 10
    put!(clk.ac[1].forth, DiscreteEvents.Run(Δt, false))
    tns[1] += take!(clk.ac[1].back).t
    iter[1] += 1
end
t1 = time_ns() - t1
println("$iter ticks took $(Int(t1)*1e-9) s, clock time $(Int(tns[1])*1e-9) s")
@test c2.clock.time ≈ 10
@test c2.clock.scount == 1000
@test a[1] == 2
@test b[1] == 2
@test c[1] == 2000

println("... parallel clock 2 reset ...")
put!(clk.ac[1].forth, DiscreteEvents.Reset(true))
@test take!(clk.ac[1].back).x == 1

println("... testing register(!) five cases ... ")
ev1 = DiscreteEvents.DiscreteEvent(fun(incr!, a), 1.0, 0.0, 1)
ev2 = DiscreteEvents.DiscreteEvent(fun(incr!, a), 2.0, 0.0, 1)
ev3 = DiscreteEvents.DiscreteEvent(fun(incr!, a), 3.0, 0.0, 1)
ev4 = DiscreteEvents.DiscreteEvent(fun(incr!, a), 4.0, 0.0, 1)
DiscreteEvents._register(clk, ev1, 1)              # 1. register ev1 to clk
@test DiscreteEvents._nextevent(clk) == ev1
DiscreteEvents._register(clk, ev1, 2)              # 2. register ev1 to 1st parallel clock
sleep(sleeptime)
@test DiscreteEvents._nextevent(c2.clock) == ev1
DiscreteEvents._register(c2, ev2, 2)               # 3. register ev2 directly to 1st parallel clock
sleep(sleeptime)
@test length(c2.clock.sc.events) == 2
DiscreteEvents._register(c2, ev2, 1)               # 4. register ev2 back to master
sleep(sleeptime)
@test length(clk.sc.events) == 2
DiscreteEvents._register(c2.clock, ev3, 1)
sleep(sleeptime)
@test length(clk.sc.events) == 3
clk.state = DiscreteEvents.Busy()  # mock busy master
DiscreteEvents._register(c2.clock, ev4, 1)
clk.state = DiscreteEvents.Idle()
@test take!(clk.ac[1].back).ev == ev4

if nthreads() > 2                           # This fails on CI (only 2 threads)
    DiscreteEvents._register(c2, ev2, 3)           # 5. register ev2 to another parallel clock
    c3 = pclock(clk, 3)
    sleep(sleeptime)
    @test DiscreteEvents._nextevent(c3.clock) == ev2
end

println("... testing API on multiple threads ...")

function clocks_ok(clk::Clock)
    ok = true
    for i in eachindex(clk.ac)
        while isready(clk.ac[i].back)
            token = take!(clk.ac[i].back)
            if token isa Error
                println("clock $i error: ", token.exc)
                ok = false
            end
        end
        if istaskfailed(clk.ac[i].ref[])
            println("clock $i error, stacktrace:")
            println(clk.ac[i].ref[])
            ok = false
        end
    end
    return ok
end

clock = Vector{Clock}()
push!(clock, clk)
for i in eachindex(clk.ac)
    put!(clk.ac[i].forth, DiscreteEvents.Query())
    push!(clock, take!(clk.ac[i].back).x.clock )
end
println("    got $(length(clock)) clocks")
@test sum(length(c.sc.events) for c in clock) ≥ 4
resetClock!(clk)                            # test resetClock! on multiple clocks
@test clk.Δt == Δt                          # must still be the same
@test sum(length(c.sc.events) for c in clock) == 0
println("    resetClock!     ok")

x = zeros(Int, nthreads())
event!(clk, ()->x[1]+=1, 1)
sleep(sleeptime); @test clocks_ok(clk)
event!(clk, ()->x[2]+=1, 1, cid=1)
sleep(sleeptime); @test clocks_ok(clk)
event!(c2,  ()->x[2]+=1, 2)
sleep(sleeptime); @test clocks_ok(clk)
event!(c2,  ()->x[1]+=1, 2, cid=2)
sleep(sleeptime); @test clocks_ok(clk)
if nthreads() > 2
    event!(c3,  ()->x[3]+=1, 1, cid=3)
    sleep(sleeptime); @test clocks_ok(clk)
end
sleep(sleeptime)
@test length(clk.sc.events) == 2
@test length(c2.clock.sc.events) == 2
nthreads() > 2 && @test length(c3.clock.sc.events) == 1
error = false
@test !error
println("    distributed events ok")
run!(clk, 3)
@test x[1] == 2
@test x[2] == 2
nthreads() > 2 && @test x[3] == 1
println("    first parallel run ok")

println("... collapse! ...")
@test_warn "only the master clock" collapse!(c2.clock)
ac = clk.ac
collapse!(clk)
@test all(x->istaskdone(x.ref[]), ac)
@test all(x->x.ch.state == :closed, ac)


text:
using TravisTest
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
@test testfun()


text:
# Needs to be installed, e.g. via:
#   julia -e 'Pkg.clone("git@github.com:loleg/TableSchema.jl.git")'
using TableSchema
import TableSchema: read, validate

include("../src/DataPackage.jl")
using DataPackage

p = Package("data/cities/datapackage.json")

println("Loaded a Data Package, resources: ", length(p.resources))

r = DataPackage.get_resource(p, "cities")

data = DataPackage.read(r)
println( "The number of cells is ", length(data[:,1]) )
println( "Sum of column 1 is ", sum([ row for row in data[1] ]) )

table = DataPackage.get_table(r)
if validate(table); println("The table is valid according to the Schema"); end


text:
# indicator of a simplex

export IndSimplex

"""
**Indicator of a simplex**

    IndSimplex(a=1.0)

Returns the indicator of the set
```math
S = \\left\\{ x : x \\geq 0, ∑_i x_i = a \\right\\}.
```
By default `a=1.0`, therefore ``S`` is the probability simplex.
"""
struct IndSimplex{T <: Union{Real, Integer}} <: ProximableFunction
  a::T
  function IndSimplex{T}(a::T) where {T <: Union{Real, Integer}}
    if a <= 0
      error("parameter a must be positive")
    else
      new(a)
    end
  end
end

is_convex(f::IndSimplex) = true
is_set(f::IndSimplex) = true

IndSimplex(a::T=1.0) where {T <: Union{Real, Integer}} = IndSimplex{T}(a)

function (f::IndSimplex)(x::AbstractArray{T}) where T <: Real
  if all(x .>= 0) && abs(sum(x)-f.a) <= 1e-14
    return zero(T)
  end
  return +Inf
end

function prox!(y::AbstractArray{T}, f::IndSimplex, x::AbstractArray{T}, gamma::Real=1.0) where T <: Real
# Implements Algorithm 1 in Condat, "Fast projection onto the simplex and the l1 ball", Mathematical Programming, 158:575–585, 2016.
# We should consider implementing the other algorithms reviewed there, and the one proposed in the paper.
  n = length(x)
  p = []
  if ndims(x) == 1
    p = sort(x, rev=true)
  else
    p = sort(x[:], rev=true)
  end
  s = 0
  for i = 1:n-1
    s = s + p[i]
    tmax = (s - f.a)/i
    if tmax >= p[i+1]
      @inbounds for j in eachindex(y)
        y[j] = x[j] < tmax ? 0.0 : x[j] - tmax
      end
      return zero(T)
    end
  end
  tmax = (s + p[n] - f.a)/n
  @inbounds for j in eachindex(y)
    y[j] = x[j] < tmax ? 0.0 : x[j] - tmax
  end
  return zero(T)
end

fun_name(f::IndSimplex) = "indicator of the probability simplex"
fun_dom(f::IndSimplex) = "AbstractArray{Real}"
fun_expr(f::IndSimplex) = "x ↦ 0 if x ⩾ 0 and sum(x) = a, +∞ otherwise"
fun_params(f::IndSimplex) = "a = $(f.a)"

function prox_naive(f::IndSimplex, x::AbstractArray{T}, gamma::Real=1.0) where T <: Real
  low = minimum(x)
  upp = maximum(x)
  v = x
  s = Inf
  for i = 1:100
    if abs(s)/f.a <= 1e-15
      break
    end
    alpha = (low+upp)/2
    v = max.(x .- alpha, 0.0)
    s = sum(v) - f.a
    if s <= 0
      upp = alpha
    else
      low = alpha
    end
  end
  return v, 0.0
end


text:
# modules prevent namespace pollution

module CoreTest

using ParserCombinator
using Test

@testset "core" begin

include("core/sources.jl")
include("core/fix.jl")
include("core/print.jl")
include("core/names.jl")
include("core/tests.jl")
include("core/slow.jl")
include("core/case.jl")
include("core/calc.jl")
include("core/debug.jl")
include("core/try.jl")

end
end

module GmlTest

using ParserCombinator
using ParserCombinator.Parsers.GML
using Test

@testset "GML" begin
include("gml/ok.jl")
include("gml/errors.jl")
include("gml/example1.jl")
include("gml/example2.jl")
# need zip files unpacking
#include("gml/celegansneural.jl")
#include("gml/polblogs.jl")
#include("gml/10k-49963.jl")

end
end

module DotTest

using ParserCombinator
using ParserCombinator.Parsers.DOT
using Test

@testset "DOT" begin
include("dot/example.jl")
include("dot/fragments.jl")
include("dot/subgraphs.jl")
include("dot/examples.jl")
end
end


text:
# TODO: Make `UniformSampling` and `Prior` algs + just use `Sampler`
# That would let us use all defaults for Sampler, combine it with other samplers etc.
"""
Robust initialization method for model parameters in Hamiltonian samplers.
"""
struct SampleFromUniform <: AbstractSampler end
struct SampleFromPrior <: AbstractSampler end

getspace(::Union{SampleFromPrior,SampleFromUniform}) = ()

# Initializations.
init(rng, dist, ::SampleFromPrior) = rand(rng, dist)
function init(rng, dist, ::SampleFromUniform)
    return istransformable(dist) ? inittrans(rng, dist) : rand(rng, dist)
end

init(rng, dist, ::SampleFromPrior, n::Int) = rand(rng, dist, n)
function init(rng, dist, ::SampleFromUniform, n::Int)
    return istransformable(dist) ? inittrans(rng, dist, n) : rand(rng, dist, n)
end

"""
    Sampler{T}

Generic sampler type for inference algorithms of type `T` in DynamicPPL.

`Sampler` should implement the AbstractMCMC interface, and in particular
[`AbstractMCMC.step`](@ref). A default implementation of the initial sampling step is
provided that supports resuming sampling from a previous state and setting initial
parameter values. It requires to overload [`loadstate`](@ref) and [`initialstep`](@ref)
for loading previous states and actually performing the initial sampling step,
respectively. Additionally, sometimes one might want to implement [`initialsampler`](@ref)
that specifies how the initial parameter values are sampled if they are not provided.
By default, values are sampled from the prior.
"""
struct Sampler{T} <: AbstractSampler
    alg::T
    selector::Selector # Can we remove it?
    # TODO: add space such that we can integrate existing external samplers in DynamicPPL
end
Sampler(alg) = Sampler(alg, Selector())
Sampler(alg, model::Model) = Sampler(alg, model, Selector())
Sampler(alg, model::Model, s::Selector) = Sampler(alg, s)

# AbstractMCMC interface for SampleFromUniform and SampleFromPrior
function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::Model,
    sampler::Union{SampleFromUniform,SampleFromPrior},
    state=nothing;
    kwargs...,
)
    vi = VarInfo()
    model(rng, vi, sampler)
    return vi, nothing
end

# initial step: general interface for resuming and
function AbstractMCMC.step(
    rng::Random.AbstractRNG, model::Model, spl::Sampler; resume_from=nothing, kwargs...
)
    if resume_from !== nothing
        state = loadstate(resume_from)
        return AbstractMCMC.step(rng, model, spl, state; kwargs...)
    end

    # Sample initial values.
    _spl = initialsampler(spl)
    vi = VarInfo(rng, model, _spl)

    # Update the parameters if provided.
    if haskey(kwargs, :init_params)
        initialize_parameters!(vi, kwargs[:init_params], spl)

        # Update joint log probability.
        # TODO: fix properly by using sampler and evaluation contexts
        # This is a quick fix for https://github.com/TuringLang/Turing.jl/issues/1588
        # and https://github.com/TuringLang/Turing.jl/issues/1563
        # to avoid that existing variables are resampled
        if _spl isa SampleFromUniform
            model(rng, vi, SampleFromPrior())
        else
            model(rng, vi, _spl)
        end
    end

    return initialstep(rng, model, spl, vi; kwargs...)
end

"""
    loadstate(data)

Load sampler state from `data`.
"""
function loadstate end

"""
    initialsampler(sampler::Sampler)

Return the sampler that is used for generating the initial parameters when sampling with
`sampler`.

By default, it returns an instance of [`SampleFromPrior`](@ref).
"""
initialsampler(spl::Sampler) = SampleFromPrior()

function initialize_parameters!(vi::AbstractVarInfo, init_params, spl::Sampler)
    @debug "Using passed-in initial variable values" init_params

    # Flatten parameters.
    init_theta = mapreduce(vcat, init_params) do x
        vec([x;])
    end

    # Get all values.
    linked = islinked(vi, spl)
    linked && invlink!(vi, spl)
    theta = vi[spl]
    length(theta) == length(init_theta) ||
        error("Provided initial value doesn't match the dimension of the model")

    # Update values that are provided.
    for i in 1:length(init_theta)
        x = init_theta[i]
        if x !== missing
            theta[i] = x
        end
    end

    # Update in `vi`.
    vi[spl] = theta
    linked && link!(vi, spl)

    return nothing
end

"""
    initialstep(rng, model, sampler, varinfo; kwargs...)

Perform the initial sampling step of the `sampler` for the `model`.

The `varinfo` contains the initial samples, which can be provided by the user or
sampled randomly.
"""
function initialstep end


text:
########### Test for LinearOperators
function test_op(A::AbstractOperator, x, y, verb::Bool = false)

  verb && (println(); show(A); println())

  Ax = A*x
  Ax2 = similar(Ax)
  verb && println("forward preallocated")
  mul!(Ax2, A, x) #verify in-place linear operator works
  verb && @time mul!(Ax2, A, x)

  @test norm(Ax .- Ax2) <= 1e-8

  Acy = A'*y
  Acy2 = similar(Acy)
  verb && println("adjoint preallocated")
  At = AdjointOperator(A)
  mul!(Acy2, At, y) #verify in-place linear operator works
  verb && @time mul!(Acy2, At, y)

  @test norm(Acy .- Acy2) <= 1e-8

  s1 = real(dot(Ax2, y))
  s2 = real(dot(x, Acy2))
  @test abs( s1 - s2 ) < 1e-8

  return Ax
end

########### Test for NonLinearOperators
function test_NLop(A::AbstractOperator, x, y, verb::Bool = false)

	verb && (println(),println(A))

	Ax = A*x
	Ax2 = similar(Ax)
	verb && println("forward preallocated")
	mul!(Ax2, A, x) #verify in-place linear operator works
	verb && @time mul!(Ax2, A, x)

	@test_throws ErrorException A'

	@test norm(Ax .- Ax2) <= 1e-8

	J = Jacobian(A,x)
	verb && println(J)

	grad = J'*y
	mul!(Ax2, A, x) #redo forward
	verb && println("adjoint jacobian mul! preallocated")
	grad2 = similar(grad)
	mul!(grad2, J', y) #verify in-place linear operator works
	verb && mul!(Ax2, A, x) #redo forward
	verb && @time mul!(grad2, J', y)

	@test norm(grad .- grad2) < 1e-8

	if all(isreal.(grad))  # currently finite difference gradient not working with complex variables
		grad3 = gradient_fd(A,Ax,x,y) #calculate gradient using finite differences

		@test norm(grad .- grad3) < 5e-4
	end

	return Ax, grad
end

############# Finite Diff for Jacobian tests

function gradient_fd(op::A,
                     y0::AbstractArray,
                     x0::AbstractArray,
                     r::AbstractArray) where {A<:AbstractOperator}

	y = copy(y0)
	J = zeros(*(size(op,1)...),*(size(op,2)...))
	h = sqrt(eps())
	for i = 1:size(J,2)
		x = copy(x0)
		x[i] = x[i]+h
		mul!(y,op,x)
		J[:,i] .= ((y.-y0)./h)[:]
	end
	return reshape(J'*r[:],size(op,2))
end

function gradient_fd(op::A,
                     y0::AbstractArray,
                     x0::ArrayPartition,
                     r::AbstractArray) where {A<:AbstractOperator}
  N = length(x0.x)
	y = copy(y0)
	grad = zero(x0)
	J =  [ zeros(*(size(op,1)...),*(sz2...)) for sz2 in size(op,2)]

	h = sqrt(eps())
	for ii in eachindex(J)
		for i = 1:size(J[ii],2)
			x = deepcopy(x0)
			x.x[ii][i] = x.x[ii][i]+h
			mul!(y,op,x)
			J[ii][:,i] .= ((y.-y0)./h)[:]
		end
		grad.x[ii] .= reshape(J[ii]'*r[:],size(op,2)[ii])
	end
	return grad
end

function gradient_fd(op::A,
                     y0::ArrayPartition,
                     x0::AbstractArray,
                     r::ArrayPartition) where {A<:AbstractOperator}
  N = length(y0.x)
	grad = zero(x0)
	y    = zero(y0)
	J = [ zeros(*(sz1...),*(size(op,2)...)) for sz1 in size(op,1)]

	h = sqrt(eps())
	for i in eachindex(x0)
		x = deepcopy(x0)
		x[i] = x[i]+h
		mul!(y,op,x)
		for ii in eachindex(J)
			J[ii][:,i] .= ((y.x[ii].-y0.x[ii])./h)[:]
		end
	end
	for ii in eachindex(J)
		grad .+= reshape(J[ii]'*r.x[ii],size(op,2))
	end
	return grad
end

function gradient_fd(op::A,
                     y0::ArrayPartition,
                     x0::ArrayPartition,
                     r::ArrayPartition) where {A<:AbstractOperator}
	grad = zero(x0)
	y    = zero(y0)
  M = length(x0.x)
  N = length(y0.x)
	J = [ zeros(*(size(op,1)[i]...),*(size(op,2)[ii]...)) for ii = 1:M, i = 1:N ]

	h = sqrt(eps())
	for i = 1:M
		for iii in eachindex(x.x[i])
			x = deepcopy(x0)
			x.x[i][iii] = x.x[i][iii]+h
			mul!(y,op,x)

			for ii = 1:N
				J[i,ii][:,iii] .= ((y.x[ii].-y0.x[ii])./h)[:]
			end
		end
	end

	for ii = 1:N, i = 1:M
		grad.x[i] .+= reshape(J[i,ii]'*r.x[ii],size(op,2)[i])
	end
  return grad
end


text:
module SimpleSparseArrays
using Base.Cartesian
using LinearAlgebra

include("core.jl")
include("testing.jl")

end#module


text:
using Images
using ImageView

img = load("proxy-image.png")
gray_img = Gray.(img) .|> Float64
inv_img = 1 .- (gray_img .> 0) .|> Bool
thin_img = thinning(inv_img)
save("thin_img.png",thin_img)
imshow(thin_img)

max_stars_repo_name: Ved-Mahajan/JavisIsAwsome

max_issues_repo_name: Ved-Mahajan/JavisIsAwsome

max_forks_repo_name: Ved-Mahajan/JavisIsAwsome


text:
using CSV, StatsPlots; pyplot()

data1 = CSV.read("../data/machine1.csv", header=false)[:,1]
data2 = CSV.read("../data/machine2.csv", header=false)[:,1]
data3 = CSV.read("../data/machine3.csv", header=false)[:,1]

boxplot([data1,data2,data3], c=[:blue :red :green], label="",
	xticks=([1:1:3;],["1", "2", "3"]), xlabel="Machine type",
	ylabel="Pipe Diameter (mm)")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook