text:
# EVs::Array{Lar.Cells,1} 	one for each 3-cell (relation Edges-Vertices)
# FVs::Array{Lar.Cells,1} 	one for each 3-cell (relation Faces-Vertices)

# julia> @show V,EVs,FVs; # triangulation of 2-faces
(V, EVs, FVs) = ([1.01181 0.215639 0.516927 0.449016 0.91979 0.123616 0.147315 0.160136 0.213681 0.156419 0.699248 0.157422 0.242039 0.325549 0.591671 0.621329 0.628539 1.02252 0.226347 0.930498 0.134324 0.195738 0.992501 0.183203 0.208451 0.513984 0.566886 0.531412 0.511538 0.41632 0.223159 0.218365 0.221725 0.225755 0.81974 0.678562 0.530265 0.508469 0.514549 0.520034 0.294128 0.440629 0.305963 0.921653 0.129033 0.129737 0.131923 0.12989 0.319954 0.159637 0.940964 0.155899 0.150091 0.159038 0.178925 0.132213 0.133682 0.209488 0.212755 0.21717 0.045831 -0.301827 0.348275 0.000617 0.165572 -0.101237 -0.178503 -0.189822 0.579367 0.23171 0.881811 0.534154 0.662595 0.59629 0.344782 0.433038 0.753477 0.716324 0.740326 0.501848 0.568661 0.459863 0.222618 0.592366 0.594961 0.207386 0.394164 0.265922 0.522271 0.530034 0.293353 -0.180789 -0.052178 0.627953 -0.190635 0.489496 0.229234 0.579478 -0.060886 0.289468 -0.079901 -0.023388 0.656742 0.126044 -0.161846 0.518285 0.526925 0.635513 0.127723 0.603362 0.496221 0.064319 0.180298 0.576293 0.525604 0.27652 -0.087513 0.52527 0.161237 0.509324 0.145291 0.758074 0.394041 0.526972 0.355908 0.399061 0.238536 0.338359 0.238295 0.27631 0.564484 0.005828 0.294002 1.01221 1.30039 0.741732 1.02991; 0.160033 0.06801 0.102833 0.094984 0.956278 0.864255 0.659194 0.548263 0.084949 0.580426 0.316997 0.583871 0.542193 0.461307 0.305687 0.514086 0.549482 0.160649 0.068627 0.956895 0.864872 0.87197 0.420396 0.441934 0.223478 0.153821 0.413502 0.239372 0.231569 0.424659 0.068443 0.068167 0.06836 0.068593 0.137692 0.121256 0.104236 0.101329 0.10214 0.102872 0.076463 0.093529 0.078046 0.956385 0.864567 0.864608 0.864734 0.864616 0.886478 0.867939 0.780928 0.585196 0.677336 0.56981 0.393343 0.835666 0.82232 0.189753 0.173182 0.121583 -0.200245 0.102199 0.417839 0.720283 0.576781 0.512131 0.056734 0.066943 -0.35354 -0.051096 0.264543 0.566987 0.45525 0.512932 0.179982 0.360343 0.259028 0.186057 0.26983 0.038311 0.374082 0.280072 0.656498 0.383734 0.4692 0.63114 0.510617 -0.094786 0.326884 0.32346 -0.106882 0.152914 0.682359 0.543901 0.059204 -0.079254 -0.026271 0.32573 0.192415 0.485962 0.144527 0.956374 0.817917 0.925953 0.333219 0.194762 0.233648 0.615857 0.837779 0.614035 -0.015242 0.10294 0.139662 0.381916 0.218367 -0.102028 0.146722 0.324834 0.573584 -0.16916 0.07959 0.257702 0.506452 0.415618 0.441015 0.337141 0.551294 0.48303 0.550965 0.452356 0.181874 1.15396 0.883481 0.816151 0.545669 1.51776 1.24728; 0.196256 0.206964 0.202911 0.203825 0.196872 0.20758 0.207421 0.207335 0.206977 0.20736 0.200311 0.207344 0.206246 0.205192 0.201757 0.2012 0.201077 0.997729 1.00844 0.998346 1.00905 1.00823 0.99793 1.00873 1.00856 1.00453 1.00362 1.00423 1.0045 1.00562 0.76985 0.411 0.662501 0.964109 0.380467 0.534691 0.38142 0.882367 0.742631 0.616551 1.0042 0.83035 0.727092 0.33635 0.612986 0.665684 0.829368 0.67719 0.811536 0.82674 0.264299 0.209682 0.56973 0.310866 0.273055 0.601118 0.595637 0.795843 0.89635 0.781236 0.067745 0.601282 -0.08555 0.447986 0.194836 0.499611 0.433829 0.451328 0.502301 1.03584 0.349005 0.882542 0.685427 0.787184 0.978526 0.933793 0.504626 0.524499 0.524582 0.735709 0.772125 0.87254 0.628801 0.258519 0.608237 0.496507 0.379068 0.151378 0.079386 0.086386 0.176108 0.651402 0.1593 0.18809 0.433316 0.462105 0.451088 0.284025 0.370218 0.228267 0.391654 0.797002 0.825792 0.803328 1.07102 1.09981 1.08271 0.355548 0.603653 0.335159 0.611076 0.661604 0.803721 0.629483 1.04042 0.1446 0.377404 0.077468 0.310272 0.580364 0.813168 0.513232 0.746036 0.661025 0.756327 0.717271 0.454961 0.390781 0.454935 0.403805 1.13971 0.7676 1.5035 0.249344 0.985248 0.613139 1.34904], Array{Array{Int64,1},1}[[[1,11,16],[1,3,11],[11,13,16],[6,7,17],[1,5,17],[5,6,17],[1,16,17],[18,23,28],[18,26,28],[19,25,26],[18,19,26],[22,24,30],[21,22,24],[18,38,41],[38,41,42],[19,34,41],[18,19,41],[1,18,38],[1,3,35],[35,38,39],[1,35,38],[6,44,45],[5,6,44],[21,22,47],[22,47,50],[5,44,51],[1,5,51],[1,18,23],[1,23,51],[19,25,59],[19,34,59],[6,7,56],[6,45,56],[21,24,47],[61,63,65],[61,65,67],[62,67,68],[62,66,68],[61,62,67],[35,39,69],[39,69,80],[41,42,69],[41,69,70],[42,69,80],[61,69,70],[61,62,70],[11,13,63],[13,63,65],[35,69,91],[61,88,91],[3,90,91],[61,88,89],[61,69,91],[3,35,91],[3,11,90],[11,89,90],[11,63,89],[61,63,89],[62,66,92],[62,70,92],[34,41,70],[34,59,92],[34,70,92],[13,16,94],[13,65,94],[67,93,95],[65,67,93],[65,93,94],[30,105,107],[24,30,105],[105,106,107],[24,47,105],[47,102,104],[47,102,105],[47,50,104],[17,93,94],[7,56,93],[7,17,93],[56,93,109],[102,104,109],[93,102,109],[67,68,95],[38,42,80],[25,26,106],[59,92,105],[25,105,106],[25,59,105],[68,93,95],[66,68,93],[66,92,105],[66,102,105],[66,93,102],[16,17,94],[38,39,80],[26,106,115],[26,28,115],[106,107,115],[88,116,118],[88,89,118],[88,91,116],[89,90,118],[91,116,118],[90,91,118],[30,107,133],[22,132,133],[107,115,131],[107,131,133],[22,30,133],[104,109,132],[22,50,104],[22,104,132],[134,136,137],[134,135,137],[28,115,131],[23,28,131],[51,134,135],[23,131,135],[23,51,135],[132,133,136],[133,136,137],[44,51,134],[44,134,136],[45,56,109],[45,132,136],[45,109,132],[44,45,136],[131,133,137],[131,135,137]],[[3,4,15],[11,14,15],[3,11,15],[9,13,14],[8,9,13],[11,13,14],[3,4,36],[35,36,39],[3,35,36],[36,39,40],[31,42,43],[31,34,41],[31,41,42],[8,9,55],[31,34,60],[34,59,60],[61,63,65],[61,65,67],[62,67,68],[62,66,68],[61,62,67],[35,39,69],[39,69,80],[41,42,69],[41,69,70],[42,69,80],[61,69,70],[61,62,70],[11,13,63],[13,63,65],[35,69,91],[61,88,91],[3,90,91],[61,88,89],[61,69,91],[3,35,91],[3,11,90],[11,89,90],[11,63,89],[61,63,89],[62,66,92],[62,70,92],[34,41,70],[34,59,92],[34,70,92],[8,13,65],[8,55,99],[8,65,99],[67,99,101],[67,97,101],[65,67,99],[67,68,97],[92,97,112],[68,92,97],[59,92,113],[59,60,113],[92,112,113],[42,43,111],[42,80,111],[66,68,92],[39,80,111],[39,40,111],[9,88,117],[9,88,89],[9,99,117],[9,55,99],[9,14,89],[31,113,121],[31,60,113],[31,120,121],[40,111,120],[36,40,120],[43,111,120],[31,43,120],[97,112,121],[97,101,117],[91,97,120],[97,120,121],[88,91,97],[88,97,117],[14,15,89],[15,89,90],[4,15,90],[4,90,91],[4,36,120],[4,91,120],[99,101,117],[112,113,121]],[[8,12,13],[8,10,12],[8,52,55],[8,10,52],[52,54,55],[53,58,59],[58,59,60],[53,54,58],[52,53,54],[66,67,68],[10,64,66],[65,66,67],[10,12,65],[10,65,66],[10,52,64],[12,13,65],[52,53,64],[59,66,92],[53,64,66],[53,59,66],[8,13,65],[8,55,99],[8,65,99],[67,99,101],[67,97,101],[65,67,99],[67,68,97],[92,97,112],[68,92,97],[59,92,113],[59,60,113],[92,112,113],[66,68,92],[54,55,99],[58,60,113],[97,101,112],[58,99,101],[54,58,99],[58,112,113],[58,101,112]],[[4,14,15],[4,9,14],[2,4,9],[2,32,37],[2,4,37],[36,37,40],[4,36,37],[31,33,43],[2,9,32],[9,32,55],[31,33,60],[88,90,91],[88,89,90],[37,96,97],[32,97,101],[32,37,97],[32,99,101],[32,55,99],[96,97,111],[33,43,97],[43,97,111],[33,60,113],[33,112,113],[33,97,112],[37,40,96],[40,96,111],[9,88,117],[9,88,89],[9,99,117],[9,55,99],[9,14,89],[31,113,121],[31,60,113],[31,120,121],[40,111,120],[36,40,120],[43,111,120],[31,43,120],[97,112,121],[97,101,117],[91,97,120],[97,120,121],[88,91,97],[88,97,117],[14,15,89],[15,89,90],[4,15,90],[4,90,91],[4,36,120],[4,91,120],[99,101,117],[112,113,121]],[[7,10,12],[7,12,13],[13,16,17],[7,13,17],[26,28,29],[25,26,29],[24,25,30],[25,29,30],[47,48,50],[24,47,48],[48,57,59],[7,56,57],[24,25,48],[53,57,59],[25,48,59],[7,10,52],[7,52,53],[7,53,57],[66,67,68],[10,64,66],[65,66,67],[10,12,65],[10,65,66],[10,52,64],[12,13,65],[52,53,64],[59,66,92],[53,64,66],[53,59,66],[13,16,94],[13,65,94],[67,93,95],[65,67,93],[65,93,94],[30,105,107],[24,30,105],[105,106,107],[24,47,105],[47,102,104],[47,102,105],[47,50,104],[17,93,94],[7,56,93],[7,17,93],[56,93,109],[102,104,109],[93,102,109],[67,68,95],[25,26,106],[59,92,105],[25,105,106],[25,59,105],[68,93,95],[66,68,93],[66,92,105],[66,102,105],[66,93,102],[16,17,94],[26,106,115],[26,28,115],[106,107,115],[48,104,109],[48,50,104],[48,57,109],[29,30,115],[30,107,115],[28,29,115],[56,57,109]],[[28,29,30],[27,28,30],[48,49,50],[46,48,49],[46,48,57],[46,56,57],[76,81,82],[72,74,76],[74,76,81],[72,74,83],[74,83,85],[83,85,87],[83,86,87],[72,76,83],[30,49,50],[27,30,103],[30,49,103],[56,108,110],[49,103,108],[46,49,56],[49,56,108],[27,28,81],[27,74,81],[74,85,108],[27,74,103],[74,103,108],[30,76,83],[29,76,82],[29,30,76],[30,50,83],[48,57,83],[57,83,86],[48,50,83],[28,29,82],[28,81,82],[85,87,110],[85,108,110],[56,57,110],[86,87,110],[57,86,110]],[[20,23,27],[20,22,27],[23,27,28],[22,27,30],[44,45,46],[44,46,49],[22,49,50],[20,22,49],[20,44,49],[20,23,51],[20,44,51],[45,46,56],[73,74,81],[73,74,85],[30,49,50],[27,30,103],[30,49,103],[56,108,110],[49,103,108],[46,49,56],[49,56,108],[27,28,81],[27,74,81],[74,85,108],[27,74,103],[74,103,108],[22,30,50],[28,51,73],[28,73,81],[23,28,51],[51,85,108],[51,73,85],[51,108,110],[45,56,110],[44,51,110],[44,45,110]],[[7,10,12],[7,12,13],[13,16,17],[7,13,17],[26,28,29],[25,26,29],[24,25,30],[25,29,30],[38,39,42],[47,48,50],[24,47,48],[48,57,59],[7,56,57],[24,25,48],[53,57,59],[25,48,59],[7,10,52],[7,52,53],[7,53,57],[10,12,52],[75,76,82],[42,75,82],[39,81,82],[39,42,82],[13,84,87],[84,85,87],[13,52,87],[52,86,87],[53,83,86],[52,53,86],[12,13,52],[59,75,76],[53,59,83],[59,76,83],[13,16,84],[24,30,50],[24,47,50],[7,56,110],[17,108,110],[7,17,110],[26,38,42],[26,42,75],[25,59,75],[25,26,75],[16,84,85],[17,85,108],[16,17,85],[26,28,38],[28,39,81],[28,38,39],[30,76,83],[29,76,82],[29,30,76],[30,50,83],[48,57,83],[57,83,86],[48,50,83],[28,29,82],[28,81,82],[85,87,110],[85,108,110],[56,57,110],[86,87,110],[57,86,110]],[[33,40,43],[32,37,40],[32,33,40],[54,55,58],[32,55,58],[32,33,60],[32,58,60],[32,37,98],[32,55,100],[32,98,100],[33,43,60],[37,40,114],[37,98,114],[54,100,119],[54,55,100],[40,114,126],[114,124,126],[40,43,126],[43,58,126],[58,125,126],[43,58,60],[100,119,127],[100,127,128],[98,124,128],[98,114,124],[98,100,128],[54,127,129],[54,119,127],[58,125,129],[54,58,129],[125,126,129],[126,129,130],[124,128,130],[124,126,130],[127,128,130],[127,129,130]],[[8,12,13],[8,10,12],[40,42,43],[39,40,42],[8,52,55],[8,10,52],[52,54,55],[53,58,59],[58,59,60],[53,54,58],[52,53,54],[10,12,52],[75,76,82],[42,75,82],[39,81,82],[39,42,82],[13,84,87],[84,85,87],[13,52,87],[52,86,87],[53,83,86],[52,53,86],[12,13,52],[59,75,76],[53,59,83],[59,76,83],[13,84,100],[8,13,100],[84,98,100],[8,55,100],[42,43,60],[59,60,75],[42,60,75],[81,85,114],[39,81,114],[84,85,114],[84,98,114],[39,40,114],[54,100,119],[54,55,100],[40,114,126],[114,124,126],[40,43,126],[43,58,126],[58,125,126],[43,58,60],[100,119,127],[100,127,128],[98,124,128],[98,114,124],[98,100,128],[54,127,129],[54,119,127],[58,125,129],[54,58,129],[76,82,126],[76,83,125],[76,125,126],[83,86,129],[83,125,129],[85,87,124],[87,124,128],[81,124,126],[81,82,126],[81,85,124],[87,127,128],[86,127,129],[86,87,127]],[[3,4,15],[11,14,15],[3,11,15],[9,13,14],[8,9,13],[11,13,14],[3,4,36],[35,36,39],[3,35,36],[36,39,40],[8,9,55],[35,71,77],[73,77,79],[35,77,78],[71,73,77],[73,78,81],[35,39,78],[39,78,81],[73,78,79],[11,73,84],[11,71,73],[11,13,84],[73,84,85],[11,35,71],[3,11,35],[13,84,100],[8,13,100],[84,98,100],[8,55,100],[81,85,114],[39,81,114],[84,85,114],[84,98,114],[39,40,114],[9,14,55],[14,55,100],[78,79,114],[36,40,114],[36,78,114],[77,79,98],[15,77,98],[79,98,114],[14,98,100],[14,15,98],[4,15,36],[15,77,78],[15,36,78],[73,81,85]],[[31,42,43],[31,34,41],[31,41,42],[31,34,60],[34,59,60],[41,42,75],[41,59,75],[34,41,59],[42,43,60],[59,60,75],[42,60,75],[31,43,60]],[[40,42,43],[39,40,42],[39,42,80],[42,43,111],[42,80,111],[39,80,111],[39,40,111],[40,43,111]],[[4,14,15],[4,9,14],[2,4,9],[2,32,37],[2,4,37],[36,37,40],[4,36,37],[2,9,32],[9,32,55],[77,78,79],[32,37,98],[32,55,100],[32,98,100],[37,40,114],[37,98,114],[9,14,55],[14,55,100],[78,79,114],[36,40,114],[36,78,114],[77,79,98],[15,77,98],[79,98,114],[14,98,100],[14,15,98],[4,15,36],[15,77,78],[15,36,78]],[[48,49,50],[46,48,49],[46,48,57],[46,56,57],[49,50,104],[46,49,104],[46,56,109],[46,104,109],[48,104,109],[48,50,104],[48,57,109],[56,57,109]],[[28,29,30],[27,28,30],[27,30,107],[27,28,107],[28,107,115],[29,30,115],[30,107,115],[28,29,115]],[[76,81,82],[72,74,76],[74,76,81],[72,74,83],[74,83,85],[83,85,87],[83,86,87],[72,76,83],[74,81,85],[123,124,126],[123,125,126],[124,127,128],[123,124,127],[123,125,129],[123,127,129],[76,82,126],[76,83,125],[76,125,126],[83,86,129],[83,125,129],[85,87,124],[87,124,128],[81,124,126],[81,82,126],[81,85,124],[87,127,128],[86,127,129],[86,87,127]],[[77,78,79],[78,79,122],[77,79,122],[77,78,122]],[[1,11,16],[1,3,11],[11,13,16],[6,7,17],[1,5,17],[5,6,17],[1,16,17],[18,23,28],[18,26,28],[19,25,26],[18,19,26],[18,38,41],[38,41,42],[19,34,41],[18,19,41],[1,18,38],[1,3,35],[35,38,39],[1,35,38],[6,44,45],[5,6,44],[5,44,51],[1,5,51],[1,18,23],[1,23,51],[19,25,59],[19,34,59],[6,7,56],[6,45,56],[35,71,77],[73,77,79],[35,77,78],[71,73,77],[73,78,81],[35,39,78],[39,78,81],[73,78,79],[41,42,75],[11,73,84],[11,71,73],[11,13,84],[73,84,85],[11,35,71],[3,11,35],[41,59,75],[34,41,59],[13,16,84],[7,56,110],[17,108,110],[7,17,110],[26,38,42],[26,42,75],[25,59,75],[25,26,75],[16,84,85],[17,85,108],[16,17,85],[26,28,38],[28,39,81],[28,38,39],[78,79,122],[77,79,122],[77,78,122],[28,51,73],[28,73,81],[23,28,51],[51,85,108],[51,73,85],[51,108,110],[45,56,110],[44,51,110],[44,45,110]],[[22,24,30],[21,22,24],[21,22,47],[22,47,50],[21,24,47],[24,30,50],[24,47,50],[22,30,50]],[[20,23,27],[20,22,27],[23,27,28],[22,27,30],[44,45,46],[44,46,49],[22,49,50],[20,22,49],[20,44,49],[20,23,51],[20,44,51],[45,46,56],[49,50,104],[27,30,107],[46,49,104],[46,56,109],[46,104,109],[27,28,107],[28,107,115],[30,107,133],[22,132,133],[107,115,131],[107,131,133],[22,30,133],[104,109,132],[22,50,104],[22,104,132],[134,136,137],[134,135,137],[28,115,131],[23,28,131],[51,134,135],[23,131,135],[23,51,135],[132,133,136],[133,136,137],[44,51,134],[44,134,136],[45,56,109],[45,132,136],[45,109,132],[44,45,136],[131,133,137],[131,135,137]],[[38,39,42],[39,42,80],[38,42,80],[38,39,80]],[[33,40,43],[32,37,40],[32,33,40],[54,55,58],[32,55,58],[32,33,60],[32,58,60],[37,96,97],[32,97,101],[32,37,97],[32,99,101],[32,55,99],[96,97,111],[33,43,97],[43,97,111],[33,60,113],[33,112,113],[33,97,112],[37,40,96],[40,96,111],[54,55,99],[58,60,113],[40,43,111],[97,101,112],[58,99,101],[54,58,99],[58,112,113],[58,101,112]],[[31,33,43],[31,33,60],[33,43,60],[31,43,60]],[[73,74,81],[73,74,85],[74,81,85],[73,81,85]],[[88,90,91],[88,89,90],[88,116,118],[88,89,118],[88,91,116],[89,90,118],[91,116,118],[90,91,118]],[[123,124,126],[123,125,126],[124,127,128],[123,124,127],[123,125,129],[123,127,129],[125,126,129],[126,129,130],[124,128,130],[124,126,130],[127,128,130],[127,129,130]]])


text:
import FiniteDiff
import ForwardDiff
import ReverseDiff
# import Zygote

using Test

grads = let
    F(x, θ, deg) = (θ[1] - x[1])^deg + θ[2] * (x[2] - x[1]^deg)^deg

    F_idx_val(ca) = F(ca[Val(:x)], ca[Val(:θ)], ca[Val(:deg)])
    F_idx_sym(ca) = F(ca[:x], ca[:θ], ca[:deg])
    F_prop(ca) = F(ca.x, ca.θ, ca.deg)

    ca = ComponentArray(x = [1, 2], θ = [1.0, 100.0], deg = 2)

    (
        truth = [-400, 200],
        finite = FiniteDiff.finite_difference_gradient(ca -> F_prop(ca), ca).x,
        forward = ForwardDiff.gradient(ca -> F_prop(ca), ca).x,
        reverse = ReverseDiff.gradient(ca -> F_prop(ca), ca).x,
        # zygote = Zygote.gradient(ca -> F(ca), ca)[1].x,
    )
end

@test grads.finite ≈ grads.truth
@test grads.forward ≈ grads.truth
@test grads.reverse ≈ grads.truth
# @test grads.zygote ≈ grads.truth

ext: jl


text:
module Example
import Random

"""
    MyPareto(alpha::Float64,x0::Float64)

Pareto distribution for sampling. This does not depend on Distributions and
the RNG seed can be set. This may not be needed for future versions of Julia.
"""
struct MyPareto
    alpha::Float64
    x0::Float64
end
MyPareto(alpha=1.0) = MyPareto(alpha, 1.0)
Random.rand(p::MyPareto) = p.x0 * Random.rand()^(-one(p.alpha) / p.alpha)

"""
    makeparetodata()

Create a reproducible sorted array of Pareto distribution samples.
"""
function makeparetodata(alpha=0.5, seed=11)
    Random.seed!(seed)
    d = MyPareto(alpha)
    return [Random.rand(d) for _=1:10^6] |> sort!
end

end # module Example


text:
@testset "Test interpolation of outgoing_translation_matrix" begin

    particle_medium = Acoustic(3; ρ=10.0, c=80.0);
    medium = Acoustic(3; ρ=1.0, c=1.0);

    R = 3.0
    r = 1.0

    separation_ratio = 1.02
    ka = 0.1
    k = ka / r
    vol_fraction = 0.3

    basis_order = 2

    ω = k * real(medium.c)

    s1 = Specie(
        particle_medium, Sphere(r);
        number_density = vol_fraction / volume(Sphere(r)),
        exclusion_distance = separation_ratio
    );
    species = [s1]
    material = Material(Sphere(R),species);

    a12 = 2.0 * minimum(outer_radius(s) * s.exclusion_distance for s in material.species)

    tol = 1e-2

    Uinter = outgoing_translation_matrix(ω, medium, material;  basis_order = basis_order, tol = tol)

    rs = LinRange(a12,2R-a12,13);
    φs = LinRange(-pi,pi,13);
    θs = LinRange(0.0,pi,13);

    r2x = radial_to_cartesian_coordinates

    error = [
        begin
            U1 = outgoing_translation_matrix(medium, basis_order, ω, r2x([r,θ,φ]))
            norm(Uinter(r2x([r,θ,φ])) - U1) / norm(U1)
        end
    for r in rs, θ in θs, φ in φs];

    @test maximum(error) < 5 * tol

    @test mean(error) < tol

end


text:
# f(x)=tanh(x)

sf=1.0

Ftanh(x::Array{Float64,2})=(sf*tanh(x),[]);
Ftanh_inplace(value,aux,x::Array{Float64,2})=copy!(value,sf*tanh(x))

Dtanh(derivativeIDX,f_c,faux_c,grad_c,grad_n,x)=axpy!(sf,grad_c.*(1.-(f_c/sf).^2),grad_n)

if PROC=="GPU" ## TODO
    function Ftanh(x::CudaArray)
        tmp=copy(x)
        tanh!(sf,x,tmp)
        return (tmp,[]) # problem here is that I don't know how to free tmp -- this will cause memory leak I think.
    end
    Ftanh_inplace(value::CudaArray,aux,x::CudaArray)=tanh!(sf,x,value)

    Dtanh(derivativeIDX,f_c,faux_c,grad_c::CudaArray,grad_n::CudaArray,x::CudaArray)=Dtanh!(sf,grad_c,f_c,grad_n)
end

Derivative[Ftanh]=Dtanh
Inplace[Ftanh]=Ftanh_inplace

import Base.tanh
tanh(A::ADnode)=ADnode(Ftanh,A)
export tanh

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
function distanceArray( array::Array{T,N} ) where {T,N}

	distances = Array{Float32,N}(undef,size(array))
	center    = (size(array,1),size(array,2),size(array,3)) ./ 2 .+ 0.5

	@inbounds begin
	for dep in 1:size(array,3)
		dd2 = ( dep - center[3] )^2
		for col in 1:size(array,2)
			dc2 = ( col - center[2] )^2
			dcd = dc2 + dd2
			for row in 1:size(array,1)
				dr = ( row - center[1] )

				distances[row,col,dep] = sqrt( muladd( dr, dr, dcd ) )
			end
		end
	end end

	return distances
end

function distoid( input::Array{<:Real,N} ) where {N}

	distances = distanceArray( input )

	# distoid = sum( distances .* ( input ./ sumI ) ) = sum( distances .* input ) / sumI
	distoid = 0.0
	sumI    = 0.0
	@inbounds @simd for idx in 1:length(input)
		sumI    += input[idx]
		distoid += distances[idx]*input[idx]
	end

	return distoid / sumI
end

function distoid( input::Array{<:Real,N}, distances ) where {N}

	# distoid = sum( distances .* ( input ./ sumI ) ) = sum( distances .* input ) / sumI
	distoid = 0.0
	sumI    = 0.0
	@inbounds @simd for idx in 1:length(input)
		sumI    += input[idx]
		distoid += distances[idx]*input[idx]
	end

	return distoid / sumI
end


text:
"""
    struct MolecularData{T}

Structure describing molecular Hamiltonian. This stores the specification
of the problem and the resulting integrals and nuclear repulsion energy.
"""
struct MolecularData{T}
    spec::MolecularSpec
    nuclear_repulsion::Float64
    one_body_integrals::Array{T, 2}
    two_body_integrals::Array{T, 4}
end

## Remember `deepcopy` uses a good fallback. This is only shallow
function Base.copy(md::MolecularData)
    return MolecularData(md.spec, md.nuclear_repulsion,
                  md.one_body_integrals,
                  md.two_body_integrals)
end

function Base.:(==)(md1::MolecularData, md2::MolecularData)
    return md1.spec == md2.spec &&
        md1.nuclear_repulsion == md2.nuclear_repulsion &&
        md1.one_body_integrals == md2.one_body_integrals &&
        md1.two_body_integrals == md2.two_body_integrals
end

function ZChop.zchop!(md::MolecularData, zeps::Real=ZChop.ZEPS)
    ZChop.zchop!(md.one_body_integrals, zeps)
    ZChop.zchop!(md.two_body_integrals, zeps)
    return md
end

function ZChop.zchop(md::MolecularData, zeps::Real=ZChop.ZEPS)
    return MolecularData(md.spec,
                         ZChop.zchop(md.nuclear_repulsion, zeps),
                         ZChop.zchop(md.one_body_integrals, zeps),
                         ZChop.zchop(md.two_body_integrals, zeps))
end


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

mutable struct DebugSetting <: SwaggerModel
    detailLevel::Any # spec type: Union{ Nothing, String } # spec name: detailLevel

    function DebugSetting(;detailLevel=nothing)
        o = new()
        validate_property(DebugSetting, Symbol("detailLevel"), detailLevel)
        setfield!(o, Symbol("detailLevel"), detailLevel)
        o
    end
end # type DebugSetting

const _property_map_DebugSetting = Dict{Symbol,Symbol}(Symbol("detailLevel")=>Symbol("detailLevel"))
const _property_types_DebugSetting = Dict{Symbol,String}(Symbol("detailLevel")=>"String")
Base.propertynames(::Type{ DebugSetting }) = collect(keys(_property_map_DebugSetting))
Swagger.property_type(::Type{ DebugSetting }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_DebugSetting[name]))}
Swagger.field_name(::Type{ DebugSetting }, property_name::Symbol) =  _property_map_DebugSetting[property_name]

function check_required(o::DebugSetting)
    true
end

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


text:
@testset "Compare to MLJLinearModels.jl"  begin

@info "$(now()) Starting test set: Compare to MLJLinearModels.jl"

################################################################################
# This code is from MLJLinearModels.jl/test/testutils.jl

sparsify!(θ, s, r) = (θ .*= (rand(r, length(θ)) .< s))

function multi_rand(Mp, r)
    # Mp[i, :] sums to 1
    n, c = size(Mp)
    be   = reshape(rand(r, length(Mp)), n, c)
    y    = zeros(Int, n)
    @inbounds for i in eachindex(y)
        rp = 1.0
        for k in 1:c-1
            if (be[i, k] < Mp[i, k] / rp)
                y[i] = k
                break
            end
            rp -= Mp[i, k]
        end
    end
    y[y .== 0] .= c
    return y
end

function generate_multiclass(n, p, c; seed=53412224, sparse=1)
    r   = StableRNG(seed)
    X   = randn(r, n, p)
    X1  = MLJLinearModels.augment_X(X, true)
    θ   = randn(r, p * c)
    θ1  = randn(r, (p+1) * c)
    sparse < 1 && begin
        sparsify!(θ, sparse, r)
        sparsify!(θ1, sparse, r)
    end
    y   = zeros(Int, n)
    y1  = zeros(Int, n)

    P   = MLJLinearModels.apply_X(X, θ, c)
    M   = exp.(P)
    Mn  = M ./ sum(M, dims=2)
    P1  = MLJLinearModels.apply_X(X, θ1, c)
    M1  = exp.(P1)
    Mn1 = M1 ./ sum(M1, dims=2)

    y  = multi_rand(Mn, r)
    y1 = multi_rand(Mn1, r)
    return ((X, y, θ), (X1, y1, θ1))
end

################################################################################
# Generate data
n = 1000
p = 5
c = 3
((X0, y0, theta0), (X1, y1, theta1)) = generate_multiclass(n, p, c)  # Data without/with intercept

# Split into train/test data
ntrain = floor(Int, 0.8 * n)
ntest  = n - ntrain
rp     = randperm(n)
itrain = rp[1:ntrain]
itest  = rp[(ntrain + 1):n]
ytrain, ytest = y1[itrain], y1[itest]
Xtrain, Xtest = X1[itrain, :], X1[itest, :]

################################################################################
# Unregularized fit

# MultinomialRegression
B  = fit(ytrain, Xtrain).params
p  = [predict(B, Xtest[i, :]) for i = 1:ntest]
LL = sum(log.([p[i][ytest[i]] for i = 1:ntest]))

# MLJLinearModels
B_mlj  = mlj_fit(ytrain, Xtrain; fit_intercept=false, penalty=:none)
p_mlj  = mlj_predict(Xtest, B_mlj)
LL_mlj = sum(log.([p_mlj[i, ytest[i]] for i = 1:ntest]))
@test isapprox(LL, LL_mlj; rtol=0.01)

################################################################################
# L1 fit

# MultinomialRegression
B  = fit(ytrain, Xtrain, L1(0.5)).params
p  = [predict(B, Xtest[i, :]) for i = 1:ntest]
LL = sum(log.([p[i][ytest[i]] for i = 1:ntest]))

# MLJLinearModels
B_mlj  = mlj_fit(ytrain, Xtrain; fit_intercept=false, gamma=0.5)
p_mlj  = mlj_predict(Xtest, B_mlj)
LL_mlj = sum(log.([p_mlj[i, ytest[i]] for i = 1:ntest]))
@test isapprox(LL, LL_mlj; rtol=0.01)

################################################################################
# L2 fit

# MultinomialRegression
B  = fit(ytrain, Xtrain, L2(0.5)).params
p  = [predict(B, Xtest[i, :]) for i = 1:ntest]
LL = sum(log.([p[i][ytest[i]] for i = 1:ntest]))

# MLJLinearModels
B_mlj  = mlj_fit(ytrain, Xtrain; fit_intercept=false, lambda=0.5)
p_mlj  = mlj_predict(Xtest, B_mlj)
LL_mlj = sum(log.([p_mlj[i, ytest[i]] for i = 1:ntest]))
@test isapprox(LL, LL_mlj; rtol=0.01)

end


text:
import KernelFunctions: kernelmatrix

export pskernel, psmkernel
export getparam
###
# simple constructor for a parametric scalar kernel
###
function pskernel(w, ker::K = SqExponentialKernel()) where K <: Kernel
    l = 1 ./ w[2:end]
    w[1]*TransformedKernel(ker,ARDTransform(l))
end

function psmkernel(dims::Int)
    w = ones(dims+1) + 0.4*(rand(dims+1).-0.5)
    return uncoupledMKernel(pskernel(w), Diagonal(ones(dims)) )
end

function psmkernel(w::Array{T,1}) where T <: Real
    dims = length(w) - 1
    return uncoupledMKernel(pskernel(w), Diagonal(ones(dims)) )
end

#####
# create new kernel with updated parameters
# currently limited functions
#####
function (ker::ScaledKernel)(w)
    typeof(ker)( ker.kernel(w[2:end]),  [w[1]])
end

# should include a check that w is as long as the current parametera
function (ker::TransformedKernel)(w)
    # not sure I need this, (eventually) just used for optimization
    l = 1 ./ w
#     l = w
    typeof(ker)(ker.kernel, typeof(ker.transform)(l))
end

# ToDO: eventually might need some getfield magic, for recursion. Problem is to know when to stop ( could probably dispatch on SimpleKernel (doesn't have parameters))
# SqExponentialKernel <: KernelFunctions.SimpleKernel > true
# but also, might not want to change the parameters in say gamma exponential kernel

#####
# Obtain scalar kernel parameters
#####
function getparam(x::Array{T, 1}) where {T <: Real}
     return x
end
function getparam(ker::TransformedKernel{K, ARDTransform{Array{Float64,1}}}) where K <: Kernel
    return 1/sqrt(2.) ./ ker.transform.v
end

function getparam(ker)
    res = getparam.(getfield.(Ref(ker), fieldnames(typeof(ker))))
    return reduce(vcat, reverse(res))
end

function getparam(ker::K) where K <: KronMatrixKernel
    kernels = ker.kernels
    params = getparam.(kernels)
    return reduce(vcat, params)
end

function getparam(mker::K) where K <: MatrixKernel
    return getparam(mker.kernel)
end

#######
# compute kernel matrices
#######

###
# for Kronecker compatible kernels

function kernelmatrix(mker::K, Z) where K<:KronMatrixKernel
    # remarkably faster than mapreduce( (k,q) -> kron(kernelmatrix(k, test), q), +, K,Q)
    KE = mker.kernels
    Q = mker.Q
    mapreduce(x -> kron(kernelmatrix(x[1], Z), x[2]), +, zip(KE,Q))
end

function kernelmatrix(mker::K, a, b) where K<:KronMatrixKernel
    # remarkably faster than mapreduce( (k,q) -> kron(kernelmatrix(k, test), q), +, K,Q)
    KE = mker.kernels
    Q = mker.Q
    mapreduce(x -> kron(kernelmatrix(x[1], a, b), x[2]), +, zip(KE,Q))
end

###
# for kernels that return a (full) matrix
function kernelmatrix(mker::K, Z) where K <:MatrixKernel
    innerN= size(Z[1], 1)
    outerN = length(Z)
    n = innerN* outerN
    B = Zygote.Buffer([Z[1][1]], n,n)
    @inbounds for j in 1:outerN
        modj = 1 + (j-1) * innerN
        for i in j:outerN
            modi = 1 + (i-1) * innerN
            B[modi:modi+innerN-1, modj:modj+innerN-1] = mker(Z[i], Z[j])
        end
        for i in 1:j-1
            modi = 1 + (i-1) * innerN
            B[modi:modi+innerN-1, modj:modj+innerN-1] = B[modj:modj+innerN-1, modi:modi+innerN-1]'
        end
    end
    return copy(B)
end

function kernelmatrix(mker::K, a, b) where K <:MatrixKernel
    innerN= size(a[1], 1)
    outerN1 = length(a)
    outerN2 = length(b)
    n1 = innerN* outerN1
    n2 = innerN* outerN2
    B = Zygote.Buffer([a[1][1]], n1,n2)
    @inbounds for j in 1:outerN2
        modj = 1 + (j-1) * innerN
        for i in 1:outerN1
            modi = 1 + (i-1) * innerN
            B[modi:modi+innerN-1, modj:modj+innerN-1] = mker(a[i], b[j])
        end
    end
    return copy(B)
end

####
## Old
###

function computeK(Z, kernelT::S) where {S <:matrixkernel }
    tmp = pairwise((z1, z2) -> kernelfunctionf(z1, z2, kernelT), Z,  Symmetric)
    K = symblockreduce(tmp)
end

function computeK(a, b, kernelT::S) where {S <:matrixkernel }
    tmp = [GaussianProcessODEs.kernelfunctionf(z1, z2, kernelT) for z1 in a, z2 in b]
    K = blockreduce(tmp)
end

function computeK(Z, kernelT::S) where {S <: scalarkernel}
    tmp = pairwise((z1, z2) -> kernelfunctionf(z1, z2, kernelT), Z,  Symmetric)
    x = Z[1]
    Id = Matrix{Float64}(LinearAlgebra.I, length(x), length(x))
    kron(tmp, Id)
end

function computeK(a, b, kernelT::S) where {S <: scalarkernel}
    tmp = [GaussianProcessODEs.kernelfunctionf(z1, z2, kernelT) for z1 in a, z2 in b]
    x = a[1]
    Id = Matrix{Float64}(LinearAlgebra.I, length(x), length(x))
    kron(tmp, Id)
end

#TODO add matrix valued version
# might have to take the symblock approach, allocate whole matrix and fill block by block

function computeK(a, kernelT)
    computeK(a, a, kernelT)
end

# TODO might be obsolete
# function Kx(x, Z, ker::T) where T <: scalarkernel
#     Id = Matrix{Float64}(LinearAlgebra.I, length(x), length(x))
#     scalars = kernelfunctionf.(Ref(x), Z, Ref(ker))
#     (kron(scalars', Id))
# end

# function Kx(x, Z, ker::T) where T <: matrixkernel
#     computeK([x], Z, ker)
# end

# function dKx(x, Z, ker::T) where T <: scalarkernel
#     Id = Matrix{Float64}(LinearAlgebra.I, length(x), length(x))

#     tmp = GaussianProcessODEs.derivativekernelfunctionf.(Ref(x), Z, Ref(ker))
#     tmp = reduce(hcat, tmp)
#     (kron(tmp, Id))
# end

###
# reduce matrix of arrays into blockmatrix
###

# export symblockreduce
# this one might/should be more efficient due to
function symblockreduce(A)
    innerN = size(A[1], 1)
    outerN = size(A, 1)
    n = innerN* outerN
    B = Symmetric(Matrix{typeof(A[1][1])}(undef, n, n))
    # B = Matrix{typeof(A[1][1])}(undef, n, n)
    # @inbounds
    for j in 1:outerN
        modj = 1 + (j-1) * innerN
        for i in 1:outerN
            modi = 1 + (i-1) * innerN
            B.data[modi:modi+innerN-1, modj:modj+innerN-1] .= A[i,j]
        end
    end
    return B
end

function blockreduce(A)
    reduce(vcat, [reduce(hcat, A[i, :]) for i in 1:size(A, 1)])
end


text:
# TODO: * generalize products:
#           * everything uses offsetframes
#           * add flat product -- ModifiedManifold that prevents offsetframes from forming
#           * add default flat offset frame
#       * AbstractManifold{D, C} -- codimension
#       * ParametricCurve --> ParametricManifold{D, C}, curves are special cases
module Manifolds

using LinearAlgebra
using Random
using ForwardDiff: Dual, partials, value
using RecipesBase
using StaticArrays

include("base.jl")
include("modifiers.jl")
include("products.jl")
include("plotting.jl")

export
    # Types
    AbstractManifold,
    # Primitives
    PointSpace, Sphere, Ball, Cube, ParametricManifold,
    # Modifiers
    reparametrized, scaled, translated, flattened,
    ×
end


text:
export Chromosome,
    process,
    find_active,
    get_genes,
    clone,
    forward_connections,
    get_output_trace

abstract type Chromosome end
abstract type Node end

mutable struct CGPNode <: Node
    connections::Array{Int64}
    f::Function
    active::Bool
    p::Float64
    output::Any
end

function CGPNode(ins::Array{Int64}, f::Function)
    CGPNode(ins, f, false, 1.0, 0.0)
end

function CGPNode(ins::Array{Int64}, f::Function, active::Bool)
    CGPNode(ins, f, active, 1.0, 0.0)
end

function CGPNode(ins::Array{Int64}, f::Function, active::Bool, p::Float64)
    CGPNode(ins, f, active, p, 0.0)
end

function snap(fc::Array{Float64}, p::Array{Float64})::Array{Int64}
    map(x->indmin(abs.(p-x)), fc)
end

function process(c::Chromosome, inps::Array)::Array{Float64}
    for i in 1:c.nin
        c.nodes[i].output = inps[i]
    end
    maxlength = maximum(map(length, inps))
    for n in c.nodes
        if n.active
            output = n.output
            if n.f == Config.f_input
                output = Config.index_in(inps, (n.p+1.0)/2.0)
            else
                p = n.p
                if ~CGP.Config.weights
                    p = 1.0
                end
                output = CGP.Config.scaled(p * n.f(c.nodes[n.connections[1]].output,
                                                     c.nodes[n.connections[2]].output,
                                                     n.p))
            end
            if length(output) == 0
                output = 0.0
            elseif length(output) > maxlength
                output = output[1:maxlength]
            end
            n.output = output
        end
    end
    outs = Array{Float64}(c.nout)
    for i in eachindex(outs)
        outs[i] = mean(c.nodes[c.outputs[i]].output)
    end
    outs
end

function recur_active!(active::BitArray, connections::Array{Int64}, ind::Int64)::Void
    if ~active[ind]
        active[ind] = true
        for i in 1:2
            recur_active!(active, connections, connections[i, ind])
        end
    end
end

function find_active(nin::Int64, outputs::Array{Int64}, connections::Array{Int64})::BitArray
    active = falses(size(connections, 2)+nin)
    active[1:nin] = true
    for i in eachindex(outputs)
        recur_active!(active, connections, outputs[i])
    end
    active[1:nin] = false
    active
end

function find_active(outputs::Array{Int64}, connections::Array{Int64})::BitArray
    active = falses(size(connections, 2))
    for i in eachindex(outputs)
        recur_active!(active, connections, outputs[i])
    end
    active
end

function reset!(c::Chromosome)
    for n in c.nodes
        n.output = 0.0
    end
end

function get_positions(c::Chromosome)
    error("Must be implemented in subclass")
end

function recur_output_trace(c::Chromosome, ind::Int64, visited::Array{Int64})
    if ~contains(==, visited, ind)
        append!(visited, [ind])
        if (ind > c.nin) && (c.nodes[ind].f != Config.f_input)
            for i in c.nodes[ind].connections
                recur_output_trace(c, i, visited)
            end
        end
    end
    visited
end

function get_output_trace(c::Chromosome, output_ind::Int64)
    # similar to decode, just return a list of node indices that determine the output
    recur_output_trace(c, c.outputs[output_ind], Array{Int64}(0))
end

function get_output_trace(c::Chromosome, outputs::Array{Int64})
    if length(outputs) > 0
        return unique(reduce(vcat, map(x->get_output_trace(c, x), outputs)))
    else
        return Array{Int64}(0)
    end
end

function get_output_trace(c::Chromosome)
    # same as indexing over active nodes
    get_output_trace(c, collect(1:c.nout))
end

function node_genes(c::Chromosome)
    # number of genes per node (position, x, y, f, param)
    error("Must be implemented in subclass")
end

function get_genes(c::Chromosome, node_id::Int64)
    c.genes[(c.nin+c.nout)+((node_id-1-c.nin)*node_genes(c))+(1:node_genes(c))]
end

function get_genes(c::Chromosome, nodes::Array{Int64})
    if length(nodes) > 0
        return reduce(vcat, map(x->get_genes(c, x), nodes))
    else
        return Array{Int64}(0)
    end
end

function set_genes!(c::Chromosome, node_id::Int64, genes::Array{Float64})
    c.genes[(c.nin+c.nout)+((node_id-1-c.nin)*node_genes(c))+(1:node_genes(c))] = genes
end

function forward_connections(c::Chromosome)
    connections = [[i] for i in 1:length(c.nodes)]
    for ci in eachindex(c.nodes)
        for i in 1:2
            conn = c.nodes[ci].connections[i]
            if conn > 0
                if ~(contains(==, connections[ci], conn))
                    append!(connections[ci], [conn])
                end
                for j in eachindex(connections[conn])
                    if ~(contains(==, connections[ci], j))
                        append!(connections[ci], [j])
                    end
                end
            end
        end
    end
    connections
end

function clone(c::Chromosome)
    # call gene constructor method
    typeof(c)(deepcopy(c.genes), c.nin, c.nout)
end


text:
# # Point Vortex Model
# This example demonstrates how to implement a custom agent stepper utilising higher-order timestep integration in
# a 'ContinuousSpace' by solving the infinite plane point vortex model. It also shows some methods to manipulate and
# extract data from the DataFrame object returned by 'Agents.run!'

# The PVM is an idealised model to study the dynamics of point vortices in a 2D incompressible inviscid fluid. Typically,
# it is used to investigate quantum vortex dynamics in 2D Bose-Einstein condensates. Each pointlike vortex has an associated
# polarity and imparts a force on all other vortices with a force that falls off quadratically with distance.

using Agents
using Random
using InteractiveDynamics
using CairoMakie
using Query

mutable struct PointVortex <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    polarity::Int
end

function move_agent_rk4!(
    agent::A,
    model::ABM{<:ContinuousSpace,A},
    dt::Real = 1.0,
) where {A<:AbstractAgent}
    model.space.update_vel!(agent, model)
    # Runge-Kutta 4th order integration scheme
    K1 = agent.vel
    K2 = agent.vel .+ 0.5.*dt.*K1
    K3 = agent.vel .+ 0.5.*dt.*K2
    K4 = agent.vel .+ dt.*K3
    pos = agent.pos .+ dt .* (K1 .+ 2.0.*K2 .+ 2.0.*K3 .+ K4)./6.0
    move_agent!(agent, pos, model)
    return agent.pos
end

# Generate alternating polarity vortices around a circle with random jitter
function NPointsCircle(model, npts, extent)
    x,y,Γ = zeros(npts),zeros(npts),zeros(npts)
    for j=0:npts-1
    x[j+1] = (model.properties["sigma"]+model.properties["η"]*randn(model.rng,1)[1])*cos(j*2π/(npts)) + extent[1]/2
    y[j+1] = (model.properties["sigma"]+model.properties["η"]*randn(model.rng,1)[1])*sin(j*2π/(npts)) + extent[2]/2
    Γ[j+1] = (-1)^j
    end
    return Tuple(collect(zip(x,y))), Γ
end

function initialize(;
    numagents = 50, # Number of point vortices
    seed = 1234, # Seed for reproducibility
    properties = Dict("dt" => 0.001, "sigma" => 10.0, "η" => 2),
    extent = (2500,2500), # Bounding box (simulation fails if vortex exits)
    spacing = min(extent...) /200, # Setup for Continuous space
    )
    rng = Random.MersenneTwister(seed)
    space = ContinuousSpace(extent, spacing)
    model = ABM(PointVortex, space; rng, properties)
    posits, circ = NPointsCircle(model,numagents,extent)
    # Assign positions and polaraties to velocities, with zero velocity.
    for n in 1:numagents
        vel = Tuple(zeros(2))
        pos = posits[n]
        polarity = circ[n]
        add_agent!(pos, model, vel, polarity)
    end
    return model
end

function agent_step!(PointVortex, model)
    # For each vortex, sum contributing forces from all other vortices (excluding self)
    c = (1/sqrt(2π)) # Quantum of circulation
    for a in allagents(model)
        a.vel = (0.0,0.0)
            for b in allagents(model)
            if a.id != b.id
                a.vel = a.vel .+ (-c*b.polarity*(a.pos[2] - b.pos[2])/edistance(a,b,model)^2,
                                   c*b.polarity*(a.pos[1] - b.pos[1])/edistance(a,b,model)^2)
            end
        end
    end
    move_agent_rk4!(PointVortex, model, model.properties["dt"])
    return
end

time = 10   # Time to propagate over
dt = 0.001  # Timestep to use in RK4 integration
σ = 3       # Initial seperation of cluster
η = 1       # Randomness amplitude to cluster

steps = Int(round(time/dt))
props = Dict("dt" => dt, "sigma" => σ, "η" => η)
model = initialize(properties=props) # Initialise model
adata = [:pos, :polarity] # Choose traits to track through simulation data
agent_df, _ = Agents.run!(model, agent_step!, steps; adata)

# Plot initial positions of each vortex as a green cross
fig,ax = scatter(agent_df.pos[1:model.agents.count];
    color=:green,
    marker=:cross,
    label="",
    show_axis = false
)

# Remove Decorations from plot for cleaner look
hidedecorations!(ax)  # hides ticks, grid and lables
hidespines!(ax)  # hide the frame

# For each vortex, query its path throughout simulation and plot it
for j=1:model.agents.count

# Harvest negative vortex paths
vpathN =  @from i in agent_df begin
            @where i.id == j
            @where i.polarity == -1
            @select i.pos
            @collect
            end

# Harvest positive vortex paths
vpathP =  @from i in agent_df begin
            @where i.id == j
            @where i.polarity == 1
            @select i.pos
            @collect
            end

# Plot negative vortex path with varying alpha to better illustrate motion
plot!(vpathN;
color=tuple.(:red, sin.(LinRange(0,π/2,steps+1))),
markersize=1,
label=""
)

# Plot positive vortex path with varying alpha to better illustrate motion
plot!(vpathP;
color=tuple.(:blue, sin.(LinRange(0,π/2,steps+1))),
markersize=1,
label=""
)

end

# Harvest final positions of negative vortices
EptsN =  @from i in agent_df begin
                        @where i.step == steps
                        @where i.polarity == -1
                        @select i.pos
                        @collect
                        end

# Harvest final positions of positive vortices
EptsP =  @from i in agent_df begin
            @where i.step == steps
            @where i.polarity == 1
            @select i.pos
            @collect
            end

# Plot final positions with circulation indications
scatter!(EptsP, color=:blue, marker=:utriangle)
scatter!(EptsN, color=:red, marker=:circle)

# Display figure
display(fig)


text:
"""
resnet_spec = Dict(
                18 => ("basic_block", [2, 2, 2, 2], [64, 64, 128, 256, 512]),
                34 => ("basic_block", [3, 4, 6, 3], [64, 64, 128, 256, 512]),
                50 => ("bottle_neck", [3, 4, 6, 3], [64, 256, 512, 1024, 2048]),
               101 => ("bottle_neck", [3, 4, 23, 3], [64, 256, 512, 1024, 2048]),
               152 => ("bottle_neck", [3, 8, 36, 3], [64, 256, 512, 1024, 2048]))

function _conv3x3(channels, stride, in_channels)
   return Conv2D(channels, (3,3), strides=stride, padding=1,
                    use_bias=false)
end

function BottleneckV2(channels, stride, downsample=false, in_channels=0)
    newchannels = floor(Int, channels/4)
    seq = Sequential(
                BatchNorm(),
                Conv2D(newchannels, (1,1), strides=1, use_bias=False),
                BatchNorm(),
                _conv3x3(newchannels, stride, newchannels),
                BatchNorm(),
                Conv2D(channels, (1,1), strides=1, use_bias=False),
            )

    if downsample
       add(seq, Conv2D(channels, (1,1), stride, use_bias=false))
    else
       self.downsample = None
    end

end

function forward(self::BottleneckV2, F, x)
       residual = x
       x = self.bn1(x)
       act = F.npx.activation if is_np_array() else F.Activation
       x = act(x, act_type='relu')
       if self.downsample:
           residual = self.downsample(x)
       x = self.conv1(x)

       x = self.bn2(x)
       x = act(x, act_type='relu')
       x = self.conv2(x)

       x = self.bn3(x)
       x = act(x, act_type='relu')
       x = self.conv3(x)

       return x + residual
end

"""


text:
o = out1[1]
na = o[:thetas_expect].*NaN
names = []
for n in o[:theta0][:names]
    for i=1:length(n[2])
        push!(names, n[1])
    end
end

# t_e = (hcat([o isa Dict ? o[:thetas_expect] : na for o in out1]...),
#        hcat([o isa Dict ? o[:thetas_expect] : na for o in out2]...))

t_e = (hcat([o[:thetas_expect] for o in out1 if o isa Dict]...),
       hcat([o[:thetas_expect] for o in out2 if o isa Dict]...))

t_med = (hcat([hcat(o[:thetas_quartiles]...)[3,:] for o in out1 if o isa Dict]...),
         hcat([hcat(o[:thetas_quartiles]...)[3,:] for o in out2 if o isa Dict]...))

t_mode = (hcat([o[:thetas_mode] for o in out1 if o isa Dict]...),
          hcat([o[:thetas_mode] for o in out2 if o isa Dict]...))

i = 1
ps = []
toplot = t_e
for jj=1:size(t_e[1],1)
    push!(ps, Plots.histogram(toplot[i][jj,:], lw=0, markershape=:circle, reuse=false, label=names[jj]))
#    display(ps[end])
end
Plots.plot(ps..., layout=(size(t_e[1],1),1), reuse=false)

# Looks ok for Svalbard mar-15


text:
# =========================================================================================
"""
    writevtk{T, M <: Model{T}}(
        stream::IOStream,
        model ::M
    )

Creates a VTK-compatible output file from a given surface or volume model. The exact file
type is determined by the given model:

| Model type    | Resulting file type  |
|---------------|----------------------|
| Surface model | VTK PolyData         |
| Volume model  | VTK UnstructuredGrid |

# Specification
<https://vtk.org/VTK/img/file-formats.pdf>

# Return type
`Void`

# Alias

    writevtk{T, M <: Model{T}}(
        fname::String,
        model::M
    )

Creates the VTK file by name rather than `IOStream` object.
"""
function writevtk(
        stream::IOStream,
        model ::Model{T, Triangle{T}}
    ) where T
    xdoc = XMLDocument()
    revidx = reverseindex(model.nodes)

    # VTKFile
    xroot = create_root(xdoc, "VTKFile")
    set_attribute(xroot, "type", "PolyData")

    # PolyData/Piece
    xpiece = new_child(new_child(xroot, "PolyData"), "Piece")
    set_attribute(xpiece, "NumberOfPoints", length(model.nodes))
    set_attribute(xpiece, "NumberOfVerts",  0)
    set_attribute(xpiece, "NumberOfLines",  0)
    set_attribute(xpiece, "NumberOfStrips", 0)
    set_attribute(xpiece, "NumberOfPolys", length(model.elements))

    # Points
    xpoints = new_child(new_child(xpiece, "Points"), "DataArray")
    set_attribute(xpoints, "type", "$T")
    set_attribute(xpoints, "NumberOfComponents", "3")
    set_attribute(xpoints, "format", "ascii")
    add_text(xpoints, join(unpack(model.nodes), " "))

    # Polys
    xpolys = new_child(xpiece, "Polys")

    # Polys/connectivity
    xconn = new_child(xpolys, "DataArray")
    set_attribute(xconn, "type", "Int32")
    set_attribute(xconn, "Name", "connectivity")
    set_attribute(xconn, "format", "ascii")
    add_text(xconn, join([revidx[objectid(n)]-1 for n in
        unpack([Vector{T}[o.v1, o.v2, o.v3] for o in model.elements])], " "))

    # Polys/offsets
    xoffs = new_child(xpolys, "DataArray")
    set_attribute(xoffs, "type", "Int32")
    set_attribute(xoffs, "Name", "offsets")
    set_attribute(xoffs, "format", "ascii")
    add_text(xoffs, join([3i for i in 1:length(model.elements)], " "))

    # create .vtp file
    println(stream, string(xdoc))
    nothing
end

function writevtk(
        stream::IOStream,
        model ::Model{T, Tetrahedron{T}}
    ) where T
    xdoc = XMLDocument()
    revidx = reverseindex(model.nodes)

    # VTKFile
    xroot = create_root(xdoc, "VTKFile")
    set_attribute(xroot, "type", "UnstructuredGrid")

    # UnstructuredGrid/Piece
    xpiece = new_child(new_child(xroot, "UnstructuredGrid"), "Piece")
    set_attribute(xpiece, "NumberOfPoints", length(model.nodes))
    set_attribute(xpiece, "NumberOfCells",  length(model.elements))

    # Points
    xpoints = new_child(new_child(xpiece, "Points"), "DataArray")
    set_attribute(xpoints, "type", "$T")
    set_attribute(xpoints, "NumberOfComponents", "3")
    set_attribute(xpoints, "format", "ascii")
    add_text(xpoints, join(unpack(model.nodes), " "))

    # Cells
    xcells = new_child(xpiece, "Cells")

    # Cells/connectivity
    xconn = new_child(xcells, "DataArray")
    set_attribute(xconn, "type", "Int32")
    set_attribute(xconn, "Name", "connectivity")
    set_attribute(xconn, "format", "ascii")
    add_text(xconn, join([revidx[objectid(n)]-1 for n in
        unpack([Vector{T}[o.v1, o.v2, o.v3, o.v4] for o in model.elements])], " "))

    # Cells/offsets
    xoffs = new_child(xcells, "DataArray")
    set_attribute(xoffs, "type", "Int32")
    set_attribute(xoffs, "Name", "offsets")
    set_attribute(xoffs, "format", "ascii")
    add_text(xoffs, join([4i for i in 1:length(model.elements)], " "))

    # Cells/types
    xtypes = new_child(xcells, "DataArray")
    set_attribute(xtypes, "type", "UInt8")
    set_attribute(xtypes, "Name", "types")
    set_attribute(xtypes, "format", "ascii")
    add_text(xtypes, "10 "^length(model.elements))

    # create .vtu file
    println(stream, string(xdoc))
    nothing
end

function writevtk(
        fname::String,
        model::M
    ) where {T, M <: Model{T}}
    open(fh -> writevtk(fh, model), fname, "w")
end


text:
module StaticModules

export StaticModule, @staticmodule, @with, @const_staticmodule

using JuliaVariables, MLStyle

struct StaticModule{Name, names, T <: Tuple}
    nt::NamedTuple{names, T}
    StaticModule{Name}(nt::NamedTuple{names, T}) where {Name, names, T} = begin
        @assert Name isa Symbol
        new{Name, names, T}(nt)
    end
end
Base.propertynames(::StaticModule{Name, names}) where {Name, names} = names
Base.getproperty(sm::StaticModule, s::Symbol) = getproperty(getfield(sm, :nt), s)
function Base.show(io::IO, sm::StaticModule{Name}) where {Name}
    nt = getfield(sm, :nt)
    n = !isempty(nt) ? maximum(s -> length(String(s)), keys(nt)) : 0
    print(io, "StaticModule $Name")
    if !get(io, :compact, false)
        print(io, " containing")
        foreach(keys(nt)) do k
            pad = mapreduce(_ -> " ", *, 1:(n - length(String(k))), init="")
            print(io, "\n  $pad$k = $(repr(nt[k]; context=:compact => true))")
        end
    end
end

function get_let_locals(ex::Expr)
    vars = (solve_from_local! ∘ simplify_ex)(ex).args[1].bounds
    map(x -> x.name, vars)
end

macro staticmodule(name, blck::Expr)
    blck = macroexpand(__module__, blck)
    @assert blck.head == :block
    locals = map(get_let_locals(blck)) do s
        :($s = $s)
    end
    push!(blck.args, :((; $(locals...) )))
    quote
        $name = let; $StaticModule{$(QuoteNode(name))}($blck) end
    end |> esc
end

macro const_staticmodule(name, blck::Expr)
    @assert blck.head == :block
    locals = map(get_let_locals(blck)) do s
        :($s = $s)
    end
    push!(blck.args, :((; $(locals...) )))
    quote
        const $name = let; $StaticModule{$(QuoteNode(name))}($blck) end
    end |> esc
end

_get_outers(_) = Symbol[]
_get_outers(x::Var) = x.is_global ? [x.name] : Symbol[]
function _get_outers(ex::Expr)
    @match ex begin
        Expr(:(=), _, rhs) => _get_outers(rhs)
        Expr(:tuple, _..., Expr(:(=), _, rhs)) => _get_outers(rhs)
        Expr(_, args...) => mapreduce(_get_outers, vcat, args)
    end
end

get_outers(ex) = (unique! ∘ _get_outers ∘ solve_from_local! ∘ simplify_ex)(ex)

function is_named_tuple(ex::Expr)
    (ex.head == :tuple) && ex.args[1] isa Expr && ((ex.args[1].head == :(=)) || ex.args[1].head == :parameters)
end

macro with(sm, blck::Expr)
    blck = macroexpand(__module__, blck)
    outers = get_outers(blck)
    if sm isa Expr && sm.head == :tuple && !is_named_tuple(sm)
        sms = sm.args
    else
        sms = [sm]
    end
    defs = map(outers) do s
        switch = foldr(sms, init=s) do M, ex
            :(($(QuoteNode(s)) ∈ $propertynames($M)) ? $M.$s : $ex)
        end
        :($s = $switch)
    end
    quote
        let $(defs...)
            $blck
        end
    end |> esc
end

end # module


text:
################################################################################
function loss( F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace; kwargs...)

	predictions = unique([ s.z for branch ∈ branches for s ∈ branch if s.bif ]; atol=2*step(targets.parameter))
	λ = length(targets.targets)-length(predictions)

	if λ≠0
		Φ = measure(F,branches,θ,targets)
		return errors(predictions,targets) - λ*log(Φ)
	else
		return errors(predictions,targets)
	end
end

################################################################################
function errors( predictions::AbstractVector{<:BorderedArray}, targets::StateSpace)
	return mean( p′-> mean( z->(z.p-p′)^2, predictions; type=:geometric ), targets.targets; type=:arithmetic )
end

realeigvals(J) = real(eigen(J).values)
function measure( F::Function, z::BorderedArray, θ::AbstractVector, targets::StateSpace )
	λ,dλ = realeigvals(∂Fu(F,z,θ)), derivative( z -> realeigvals(∂Fu(F,z,θ)), z, tangent_field(F,z,θ) )
	return window_function(targets.parameter,z) * mapreduce( (λ,dλ) -> 2 / ( 2 + abs(sinh(2λ)/dλ) ), +, λ, dλ )
end

# Giles, M. 2008. An extended collection of matrix derivative results for forward and reverse mode automatic differentiation
# credit: @mateuszbaran https://github.com/JuliaManifolds/Manifolds.jl/pull/27#issuecomment-521693995
import LinearAlgebra: eigen

function make_eigen_dual(val::Real, partial)
    Dual{tagtype(partial)}(val, partial.partials)
end

function make_eigen_dual(val::Complex, partial::Complex)
    Complex(Dual{tagtype(real(partial))}(real(val), real(partial).partials),
        Dual{tagtype(imag(partial))}(imag(val), imag(partial).partials))
end

function eigen(A::StridedMatrix{Dual{T,V,N}}) where {T,V,N}

    Av = map(a -> a.value, A)
    λ,U = eigen(Av)

    dΛ = U \ A * U
    dλ = diag(dΛ)

    F = similar(Av, eltype(λ))
    for i ∈ axes(A,1), j ∈ axes(A,2)
        if i == j
            F[i, j] = 0
        else
            F[i, j] = inv(λ[j] - λ[i])
        end
    end
    dU = U * (F .* dΛ)

    for i ∈ eachindex(dU)
        dU[i] = make_eigen_dual(U[i], dU[i])
    end
    return Eigen(dλ, dU)
end

################################################################################
function measure( F::Function, branch::Branch, θ::AbstractVector, targets::StateSpace )
	return sum( s -> measure(F,s.z,θ,targets)*s.ds, branch )
end

function measure( F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace )
	return sum( branch -> measure(F,branch,θ,targets), branches )
end

########################################################## bifurcation distance and velocity
distance(F::Function,z::BorderedArray,θ::AbstractVector) = [ F(z,θ); det(∂Fu(F,z,θ)) ]
function velocity( F::Function, z::BorderedArray, θ::AbstractVector; newtonOptions=NewtonPar(verbose=false,maxIter=800,tol=1e-6) )
	∂implicit, _, _ = newtonOptions.linsolver( -jacobian( z->distance(F,z,θ), z )', [zero(z.u);one(z.p)] )
	return gradient( θ->distance(F,z,θ)'∂implicit, θ )
end

################################################################################
function tangent_field( F::Function, z::BorderedArray, θ::AbstractVector )
	field = kernel(∂Fz(F,z,θ);nullity=length(z.p))
	return norm(field)^(-1) * BorderedArray( field[Not(end)], field[end] ) # unit tangent field
end

################################################################################
using SpecialFunctions: erf
function window_function( parameter::AbstractVector, z::BorderedArray; β::Real=10 )
	pMin,pMax = extrema(parameter)
	return ( 1 + erf((β*(z.p-pMin)-3)/√2) )/2 * ( 1 - ( 1 + erf((β*(z.p-pMax)+3)/√2) )/2 )
end


text:
function modalparam(A,C,dt)
	#
	# Compute modal parameters from A, C and dt
	#
	# javier.cara@upm.es, 2016-02
	#

	no,ns = size(C)

	# allocating
	wm1 = zeros(1,ns)
	zm1 = zeros(1,ns)
	vm1 = zeros(no,ns)

	# eigenvalues and eigenvectors
	D,V = eig(A)

	# C*V
	vm0 = C*V

	i = 1
	ii = 1
	fm1 = zeros(Number,no,ns) # can be real or complex numbers
	while i < ns
		d = D[i]
		if abs(imag(d)) > 1e-10 # complex number
			wm1[ii] = abs(log(d))/dt
			zm1[ii] = -real(log(d))/(wm1[ii]*dt)

			# modal shapes normalized with max component = 1
			maxval = vm0[1,i]
			for j in 2:no
				if abs(vm0[j,i]) > abs(maxval)
					maxval = vm0[j,i]
				end
			end
			vm1[:,ii] = real( vm0[:,i]/maxval )

			i = i+2
			ii = ii+1

		else
			# real eigenvalue
			i = i+1
		end
	end

	# deleting zero values
	wm = wm1[1:ii-1]
	zm = zm1[1:ii-1]
	vm = vm1[:,1:ii-1]

	# sorting frequencies
	pos=sortperm(wm)
	wm = wm[pos]
	zm = zm[pos]
	vm = vm[:,pos]

	return wm,zm,vm

end


text:
export ConvertInterval, DGP_FEOL_2

function ConvertInterval(ystar,thresholds)
    # This function translates a vector of latent variables to a vector of interval-censored outcomes.

    y = similar(ystar)
    for i in 1:length(ystar)
        y[i] = sum(ystar[i].>thresholds) + 1
    end
    return y
end
#df[:y] = ConvertInterval(df[:ystar],thresholds)

function DGP_FEOL_2(n, beta0k,thresholds1,thresholds2)

    T = 2 #hard-code T=2?
    K = 2 # also hard-coded
    nT = n*T

    beta0 = beta0k + zeros(K,1)
    D_X = Normal(0,1)
    D_u = Logistic()

    # Generate data frame.
    df = DataFrame(
        i = repeat(1:n, inner = T),
        t = repeat(1:T, outer = n),
        X1 = rand(D_X,n*T),
        X2 = rand(D_X,n*T),
        u = rand(D_u,n*T)
    );

    df_a = by(df, [:i], y-> 0.2*(mean(y[:X1]) + mean(y[:X2])))
    df_a[:a] = df_a[:x1]
    delete!(df_a,:x1)
    df = join(df,df_a,on = :i)
    df[:ystar] = 0 + df[:X1]*beta0[1] + df[:X2]*beta0[2] - df[:u] + df[:a]

    # This only works because the data is perfectly shaped and T=2.
    df[:y] = 0
    df[filter(isodd,1:nT),:y] = ConvertInterval(df[filter(isodd,1:nT),:ystar],thresholds1)
    df[filter(iseven,1:nT),:y] = ConvertInterval(df[filter(iseven,1:nT),:ystar],thresholds2)

    return df
end


text:
using MagmaThermoKinematics.Diffusion3D
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
using Plots
using LinearAlgebra
using SpecialFunctions
using Test

const CreatePlots = false      # easy way to deactivate plotting throughout

# Initialize for multiple threads (GPU is not tested here)
@init_parallel_stencil(Threads, Float64, 3);    # initialize parallel stencil in 3D

function Diffusion_Gaussian3D(Setup="3D")
    # Test the

    # Model parameters
    W,L,H                   =   300., 300., 300.;                      # Width, Length,    Height in km
    k_rock1                 =   3;
    Nx, Ny, Nz              =   100, 100, 100;                      # resolution
    Tbot                    =   0;
    SecYear                 =   365.25*24*3600
    σ                       =   15e3;                                           # halfwidth of gaussian
    Tmax                    =   1000;                                           # max. of gaussian
    TotalTime               =   3e6*SecYear;                                    # thermal cooling age
    ρ                       =   2800;                                           # Density
    cp                      =   1050;                                           # Heat capacity
    La                      =   350e3;                                          # Latent heat J/kg/K
    dx,dy,dz                =   W/(Nx-1)*1e3, L/(Ny-1)*1e3, H*1e3/(Nz-1);       # grid size [m]
    κ                       =   k_rock1./(ρ*cp);                                # thermal diffusivity
    dt                      =   min(dx^2,dy^2,dz^2)./κ/2;                       # stable timestep (required for explicit FD)

    numTime                 =   ceil(TotalTime/dt);
    dt                      =   TotalTime/numTime/20;
    nt                      =   Int(numTime);

    # Array initializations (1 - main arrays on which we can initialize properties)
    T                       =   @ones(Nx,Ny,Nz)*Tbot;
    K                       =   @ones(Nx,Ny,Nz)*k_rock1;
    Rho                     =   @ones(Nx,Ny,Nz)*ρ;
    Cp                      =   @ones(Nx,Ny,Nz)*cp;
    dPhi_dt                 =   @zeros(Nx,Ny,Nz);

    # Work array initialization
    Tnew, qx,qy,qz          =   @zeros(Nx,Ny,Nz),       @zeros(Nx-1,Ny,Nz),     @zeros(Nx, Ny-1, Nz),   @zeros(Nx,Ny,Nz-1)  # thermal solver
    Kx, Ky, Kz              =                           @zeros(Nx-1, Ny, Nz),   @zeros(Nx,Ny-1,Nz),     @zeros(Nx,Ny,Nz-1)  # thermal conductivities
    X,Y,Z                   =   @zeros(Nx,Ny,Nz),       @zeros(Nx,Ny,Nz),       @zeros(Nx,Ny,Nz)                            # 3D gridpoints

    # Set up model geometry & initial T structure
    x,y,z                   =   -W/2*1e3:dx:(-W/2*1e3+(Nx-1)*dx), -L/2*1e3:dy:(-L/2*1e3+(Ny-1)*dy), -H/2*1e3:dz:(-H/2*1e3+(Nz-1)*dz);
    coords                  =   collect(Iterators.product(x,y,z))                               # generate coordinates from 1D coordinate vectors
    X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);     # transfer coords to 3D arrays
    Grid, Spacing           =   (X,Y,Z), (dx,dy,dz);
    T                      .=   Tmax.*exp.(  -((X.^2 .+ Y.^2 .+ Z.^2)./(σ^2)) );                 # initial gaussian profile
    Tnew                   .=   T;

    #ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])
    #println("Animation directory: $(anim.dir)")

    time,time_kyrs          = 0.0, 0.0;
    err = 100;
    it = 0;
    nt = 500
    for it=1:nt

        # Perform a diffusion step
        if Setup=="3D"
            @parallel diffusion3D_step_varK!(Tnew, T, qx, qy, qz, K, Kx, Ky, Kz, Rho, Cp, dt, dx, dy, dz,  La, dPhi_dt);
            #@parallel diffusion3D_step!(Tnew, T, K, 1.0/(ρ*cp), dt, dx, dy, dz)
        end

        # diffusion in z-direction
        @parallel (1:size(T,2), 1:size(T,3)) bc3D_x!(Tnew);                                         # set lateral boundary conditions (flux-free)
        @parallel (1:size(T,1), 1:size(T,3)) bc3D_y!(Tnew);                                         # set lateral boundary conditions (flux-free)

        Tnew[:,:,1] .= Tbot; Tnew[:,:,end] .= 0.0;                                                  # bottom & top temperature (constant)

        T, Tnew         =   Tnew, T;                                                                # Update temperature
        time, time_kyrs =   time + dt, time/SecYear/1e3;                                            # Keep track of evolved time

        if mod(it,100)==0  # print progress
           # println(" Timestep $it = $(round(time/SecYear)/1e3) kyrs")
        #    x_km, z_km  =   x./1e3, z./1e3;
        #    #p1          =   heatmap(x_km, z_km, T[:,Int(Ny/2),:]',         aspect_ratio=1, xlims=(x_km[1],x_km[end]), ylims=(z_km[1],z_km[end]),   c=:inferno, title="Temperature, $(round(time_kyrs, digits=2)) kyrs",  dpi=150)
        #    p1          =   heatmap(x_km, z_km, T[Int(Nx/2),:,:]',         aspect_ratio=1, xlims=(x_km[1],x_km[end]), ylims=(z_km[1],z_km[end]),   c=:inferno, title="Temperature, $(round(time_kyrs, digits=2)) kyrs",  dpi=150)
        #
        #    plot(p1); frame(anim)
        end

    end

    x_km, z_km  =   x./1e3, z./1e3;

    # compute analytical solution

    if Setup=="3D"
        Tanal  =  Tmax./(1 + 4*time*κ/σ^2)^(3/2).*exp.(  -(X.^2 .+ Y.^2 .+ Z.^2)./(σ^2 + 4*time*κ));                     # initial gaussian profile
        fname = "Diffusion_3D_Gaussian";
    end
    Terror =  T - Tanal

    Tslice  = T[:,Int(Ny/2),:];
    Tanal1  = Tanal[:,Int(Ny/2),:];
    Terror1 = Terror[:,Int(Ny/2),:];

    if CreatePlots
        # create plot
        p1          =   heatmap(x_km, z_km, Tslice',         aspect_ratio=1, xlims=(x_km[1],x_km[end]), ylims=(z_km[1],z_km[end]),   c=:inferno, title="T  3D $(round(time_kyrs/1e3, digits=2)) Myrs",  dpi=150)
        p2          =   heatmap(x_km, z_km, Tanal1',         aspect_ratio=1, xlims=(x_km[1],x_km[end]), ylims=(z_km[1],z_km[end]),   c=:inferno, title="T anal 3D $(round(time_kyrs/1e3, digits=2)) Myrs",  dpi=150)
        p3          =   heatmap(x_km, z_km, Terror1',         aspect_ratio=1, xlims=(x_km[1],x_km[end]), ylims=(z_km[1],z_km[end]),   c=:inferno, title="T error 3D $(round(time_kyrs/1e3, digits=2)) Myrs",  dpi=150)
        plot(p1,p2,p3);
        png(fname)
    end

    error = norm(Terror[:],2)/length(Terror[:]);

    return error;        # return error
end # end of gaussian diffusion test

# Create a range of diffusion tests which calls the routines above
@testset "3D Gaussian diffusion" begin
    @test Diffusion_Gaussian3D("3D")           ≈  1.74e-5 atol=1e-4;
end;


text:
function __init__()
    R"options(width=100)"
    R"suppressWarnings(library('rstan'))"
    R"rstan_options(auto_write = TRUE)"
    R"options(mc.cores = 1)"
end


text:
import Adapt

using StaticArrays

dummy() = return

@testset "@oneapi" begin

@test_throws UndefVarError @oneapi undefined()
@test_throws MethodError @oneapi dummy(1)

@testset "low-level interface" begin
    k = zefunction(dummy)
    k()
    k(; items=1)
end

@testset "launch configuration" begin
    @oneapi dummy()

    @oneapi items=1 dummy()
    @oneapi items=(1,1) dummy()
    @oneapi items=(1,1,1) dummy()

    @oneapi groups=1 dummy()
    @oneapi groups=(1,1) dummy()
    @oneapi groups=(1,1,1) dummy()
end

@testset "launch=false" begin
    k = @oneapi launch=false dummy()
    k()
    k(; items=1)
end

@testset "inference" begin
    foo() = @oneapi dummy()
    @inferred foo()

    # with arguments, we call kernel_convert
    kernel(a) = return
    bar(a) = @oneapi kernel(a)
    @inferred bar(oneArray([1]))
end

@testset "reflection" begin
    oneAPI.code_lowered(dummy, Tuple{})
    oneAPI.code_typed(dummy, Tuple{})
    oneAPI.code_warntype(devnull, dummy, Tuple{})
    oneAPI.code_llvm(devnull, dummy, Tuple{})
    oneAPI.code_spirv(devnull, dummy, Tuple{})

    @device_code_lowered @oneapi dummy()
    @device_code_typed @oneapi dummy()
    @device_code_warntype io=devnull @oneapi dummy()
    @device_code_llvm io=devnull @oneapi dummy()
    @device_code_spirv io=devnull @oneapi dummy()

    mktempdir() do dir
        @device_code dir=dir @oneapi dummy()
    end

    @test_throws ErrorException @device_code_lowered nothing

    # make sure kernel name aliases are preserved in the generated code
    @test occursin("julia_dummy", sprint(io->(@device_code_llvm io=io optimize=false @oneapi dummy())))
    @test occursin("julia_dummy", sprint(io->(@device_code_llvm io=io @oneapi dummy())))
    @test occursin("julia_dummy", sprint(io->(@device_code_spirv io=io @oneapi dummy())))

    # make sure invalid kernels can be partially reflected upon
    let
        invalid_kernel() = throw()
        @test_throws oneAPI.KernelError @oneapi invalid_kernel()
        @test_throws oneAPI.KernelError @grab_output @device_code_warntype @oneapi invalid_kernel()
        out, err = @grab_output begin
            try
                @device_code_warntype @oneapi invalid_kernel()
            catch
            end
        end
        @test occursin("Body::Union{}", err)
    end

    # set name of kernel
    @test occursin("julia_mykernel", sprint(io->(@device_code_llvm io=io begin
        k = zefunction(dummy, name="mykernel")
        k()
    end)))
end

@testset "external kernels" begin
    @eval module KernelModule
        export external_dummy
        external_dummy() = return
    end
    import ...KernelModule
    @oneapi KernelModule.external_dummy()
    @eval begin
        using ...KernelModule
        @oneapi external_dummy()
    end

    @eval module WrapperModule
        using oneAPI
        @eval dummy() = return
        wrapper() = @oneapi dummy()
    end
    WrapperModule.wrapper()
end

@testset "calling device function" begin
    @noinline child(i) = sink(i)
    function parent()
        child(1)
        return
    end

    @oneapi parent()
end

@testset "varargs" begin
    function kernel(args...)
        oneAPI.@print(args[2])
        return
    end

    _, out = @grab_output begin
        @oneapi kernel(1, 2, 3)
        synchronize()
    end
    @test out == "2"
end

end

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

@testset "argument passing" begin

dims = (16, 16)
len = prod(dims)

@testset "manually allocated" begin
    function kernel(input, output)
        i = get_global_id(0)

        val = input[i]
        output[i] = val

        return
    end

    input = round.(rand(Float32, dims) * 100)
    output = similar(input)

    input_dev = oneArray(input)
    output_dev = oneArray(output)

    @oneapi items=len kernel(input_dev, output_dev)
    @test input ≈ Array(output_dev)
end

@testset "scalar through single-value array" begin
    function kernel(a, x)
        i = get_global_id(0)
        max = get_global_size(0)
        if i == max
            _val = a[i]
            x[] = _val
        end
        return
    end

    arr = round.(rand(Float32, dims) * 100)
    val = [0f0]

    arr_dev = oneArray(arr)
    val_dev = oneArray(val)

    @oneapi items=len kernel(arr_dev, val_dev)
    @test arr[dims...] ≈ Array(val_dev)[1]
end

@testset "scalar through single-value array, using device function" begin
    @noinline child(a, i) = a[i]
    function parent(a, x)
        i = get_global_id(0)
        max = get_global_size(0)
        if i == max
            _val = child(a, i)
            x[] = _val
        end
        return
    end

    arr = round.(rand(Float32, dims) * 100)
    val = [0f0]

    arr_dev = oneArray(arr)
    val_dev = oneArray(val)

    @oneapi items=len parent(arr_dev, val_dev)
    @test arr[dims...] ≈ Array(val_dev)[1]
end

@testset "tuples" begin
    # issue #7: tuples not passed by pointer

    function kernel(keeps, out)
        if keeps[1]
            out[] = 1
        else
            out[] = 2
        end
        return
    end

    keeps = (true,)
    d_out = oneArray(zeros(Int))

    @oneapi kernel(keeps, d_out)
    @test Array(d_out)[] == 1
end

@testset "ghost function parameters" begin
    # bug: ghost type function parameters are elided by the compiler

    len = 60
    a = rand(Float32, len)
    b = rand(Float32, len)
    c = similar(a)

    d_a = oneArray(a)
    d_b = oneArray(b)
    d_c = oneArray(c)

    @eval struct ExecGhost end

    function kernel(ghost, a, b, c)
        i = get_global_id(0)
        c[i] = a[i] + b[i]
        return
    end
    @oneapi items=len kernel(ExecGhost(), d_a, d_b, d_c)
    @test a+b == Array(d_c)

    # bug: ghost type function parameters confused aggregate type rewriting

    function kernel(ghost, out, aggregate)
        i = get_global_id(0)
        out[i] = aggregate[1]
        return
    end
    @oneapi items=len kernel(ExecGhost(), d_c, (42,))

    @test all(val->val==42, Array(d_c))
end

@testset "immutables" begin
    # issue #15: immutables not passed by pointer

    function kernel(ptr, b)
        ptr[] = imag(b)
        return
    end

    arr = oneArray(zeros(Float32))
    x = ComplexF32(2,2)

    @oneapi kernel(arr, x)
    @test Array(arr)[] == imag(x)
end

@testset "automatic recompilation" begin
    arr = oneArray(zeros(Int))

    function kernel(ptr)
        ptr[] = 1
        return
    end

    @oneapi kernel(arr)
    @test Array(arr)[] == 1

    function kernel(ptr)
        ptr[] = 2
        return
    end

    @oneapi kernel(arr)
    @test Array(arr)[] == 2
end

@testset "automatic recompilation (bis)" begin
    arr = oneArray(zeros(Int))

    @eval doit(ptr) = ptr[] = 1

    function kernel(ptr)
        doit(ptr)
        return
    end

    @oneapi kernel(arr)
    @test Array(arr)[] == 1

    @eval doit(ptr) = ptr[] = 2

    @oneapi kernel(arr)
    @test Array(arr)[] == 2
end

@testset "non-isbits arguments" begin
    function kernel1(T, i)
        sink(i)
        return
    end
    @oneapi kernel1(Int, 1)

    function kernel2(T, i)
        sink(unsafe_trunc(T,i))
        return
    end
    @oneapi kernel2(Int, 1f0)
end

@testset "splatting" begin
    function kernel(out, a, b)
        out[] = a+b
        return
    end

    out = [0]
    out_dev = oneArray(out)

    @oneapi kernel(out_dev, 1, 2)
    @test Array(out_dev)[1] == 3

    all_splat = (out_dev, 3, 4)
    @oneapi kernel(all_splat...)
    @test Array(out_dev)[1] == 7

    partial_splat = (5, 6)
    @oneapi kernel(out_dev, partial_splat...)
    @test Array(out_dev)[1] == 11
end

@testset "object invoke" begin
    # this mimics what is generated by closure conversion

    @eval struct KernelObject{T} <: Function
        val::T
    end

    function (self::KernelObject)(a)
        a[] = self.val
        return
    end

    function outer(a, val)
       inner = KernelObject(val)
       @oneapi inner(a)
    end

    a = [1f0]
    a_dev = oneArray(a)

    outer(a_dev, 2f0)

    @test Array(a_dev) ≈ [2f0]
end

@testset "closures" begin
    function outer(a_dev, val)
       function inner(a)
            # captures `val`
            a[] = val
            return
       end
       @oneapi inner(a_dev)
    end

    a = [1f0]
    a_dev = oneArray(a)

    outer(a_dev, 2f0)

    @test Array(a_dev) ≈ [2f0]
end

@testset "conversions" begin
    @eval struct Host   end
    @eval struct Device end

    Adapt.adapt_storage(::oneAPI.KernelAdaptor, a::Host) = Device()

    Base.convert(::Type{Int}, ::Host)   = 1
    Base.convert(::Type{Int}, ::Device) = 2

    out = [0]

    # convert arguments
    out_dev = oneArray(out)
    let arg = Host()
        @test Array(out_dev) ≈ [0]
        function kernel(arg, out)
            out[] = convert(Int, arg)
            return
        end
        @oneapi kernel(arg, out_dev)
        @test Array(out_dev) ≈ [2]
    end

    # convert tuples
    out_dev = oneArray(out)
    let arg = (Host(),)
        @test Array(out_dev) ≈ [0]
        function kernel(arg, out)
            out[] = convert(Int, arg[1])
            return
        end
        @oneapi kernel(arg, out_dev)
        @test Array(out_dev) ≈ [2]
    end

    # convert named tuples
    out_dev = oneArray(out)
    let arg = (a=Host(),)
        @test Array(out_dev) ≈ [0]
        function kernel(arg, out)
            out[] = convert(Int, arg.a)
            return
        end
        @oneapi kernel(arg, out_dev)
        @test Array(out_dev) ≈ [2]
    end

    # don't convert structs
    out_dev = oneArray(out)
    @eval struct Nested
        a::Host
    end
    let arg = Nested(Host())
        @test Array(out_dev) ≈ [0]
        function kernel(arg, out)
            out[] = convert(Int, arg.a)
            return
        end
        @oneapi kernel(arg, out_dev)
        @test Array(out_dev) ≈ [1]
    end
end

@testset "argument count" begin
    val = [0]
    val_dev = oneArray(val)
    for i in (1, 10, 20, 34)
        variables = ('a':'z'..., 'A':'Z'...)
        params = [Symbol(variables[j]) for j in 1:i]
        # generate a kernel
        body = quote
            function kernel(arr, $(params...))
                arr[] = $(Expr(:call, :+, params...))
                return
            end
        end
        eval(body)
        args = [j for j in 1:i]
        call = Expr(:call, :kernel, val_dev, args...)
        cudacall = :(@oneapi $call)
        eval(cudacall)
        @test Array(val_dev)[1] == sum(args)
    end
end

@testset "keyword arguments" begin
    @eval inner_kwargf(foobar;foo=1, bar=2) = nothing

    @oneapi (()->inner_kwargf(42;foo=1,bar=2))()

    @oneapi (()->inner_kwargf(42))()

    @oneapi (()->inner_kwargf(42;foo=1))()

    @oneapi (()->inner_kwargf(42;bar=2))()

    @oneapi (()->inner_kwargf(42;bar=2,foo=1))()
end

@testset "captured values" begin
    function f(capture::T) where {T}
        function kernel(ptr)
            ptr[] = capture
            return
        end

        arr = oneArray(zeros(T))
        @oneapi kernel(arr)

        return Array(arr)[1]
    end

    using Test
    @test f(1) == 1
    @test f(2) == 2
end

end

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

@testset "#55: invalid integers created by alloc_opt" begin
    function f(a)
        x = SVector(0f0, 0f0)
        v = MVector{3, Float32}(undef)
        for (i,_) in enumerate(x)
            v[i] = 1f0
        end
        a[1] = v[1]
        return nothing
    end
    @oneapi f(oneArray(zeros(Float32, 1)))
end

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


text:
include("src/core.jl")
using JLD
function createWilshire5000Data(rank, theK)
    wilshire5000Data=load("data/Wilshire5000Rank"*string(rank)*".jld", "Wilshire5000Rank"*string(rank))

    n=size(wilshire5000Data.X,2)
    μ=wilshire5000Data.μ*2769/(11*12)
    X=wilshire5000Data.X*sqrt((2769/(11*12)))
    Σ= X'*X

    min_investment=zeros(n)

    γ=(1.0/sqrt(n))*ones(n)

    # Compute r_min and r_max, in order to add minimum return constraint.

    # r_min
    m=Model(solver=MosekSolver(MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0));
    @variable(m, x[1:n]>=0.0);
    @constraint(m, sum(x)==1.0);
    @objective(m, Min,  ((1.0/(2.0*γ[1]))*x'*x +0.5*x'*Σ*x)');
    solve(m)
    r_min=μ'*getvalue(x)
    # r_max
    m2=Model(solver=MosekSolver(MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0))
    @variable(m2, x2[1:n]>=0.0)
    @constraint(m2, sum(x2)==1.0)
    @objective(m2, Max,  (μ'*x2-(1.0/(2.0*γ[1]))*x2'*x2)')
    solve(m2)
    r_max=μ'*getvalue(x2)
    r_bar=r_min+0.3*(r_max-r_min)

    A=μ'
    l=[r_bar]
    u=[1e12]

    μ=zeros(n)
    #########################################################################

    k=theK

    # A=zeros(1,n)
    # l=[0.0]
    # u=[0.0]
    f=size(X, 1)
    m=size(A, 1)
    Y=(X*X')\(X*μ)
    d=X'*Y-μ

    return SparsePortfolioData(μ, Y, d, X, A, l, u, k, n, m, f, min_investment, γ)
end

for theK in [10, 50, 100, 200]
        for rank in [100, 200, 500, 1000]
        thePortfolioData=createWilshire5000Data(rank, theK)

        sparse_Portfolio=cutting_planes_portfolios(thePortfolioData; useSOCPLB=false, minReturnConstraint=true, useHeuristic=false, useWarmStart=true, useCopyOfVariables=false, usingKelleyPrimal=true, maxCutCallbacks=0);
        @show @test(abs.(sum(sparse_Portfolio.w)-1.0)<1e-4)
        println("Variance is:", sparse_Portfolio.portfolioVariance)
        println("Return is:", sparse_Portfolio.expectedReturn)
        println("Reg term is:", 0.5*sum(sparse_Portfolio.w[i]^2/thePortfolioData.γ[sparse_Portfolio.indices[i]] for i=1:size(sparse_Portfolio.indices,1)))
        open("results/Table9ResultsOutput_withminreturn.csv","a") do fp
          println(fp, "Instance, cardinality, regularizer, method, solvetime, boundgap, boundgapSOCP, nodecount, number of iterations")
          println(fp, "rank"*string(rank), ",", thePortfolioData.k, ",", "1/sqrt(n)", ",", "OA", ",",   sparse_Portfolio.solveTime, ",", sparse_Portfolio.boundGap, ",", sparse_Portfolio.boundGapSOCP, ",", sparse_Portfolio.nodeCount, ",", sparse_Portfolio.numIters)
        end

        sparse_Portfolio_Cplex=cplex_raw_MISOCP(thePortfolioData, 600.0, 1e-4, 1)
        #@show @test(abs.(sum(sparse_Portfolio_Cplex.w)-1.0)<1e-4)
        println("Variance is:", sparse_Portfolio_Cplex.portfolioVariance)
        println("Return is:", sparse_Portfolio_Cplex.expectedReturn)
        println("Reg term is:", 0.5*sum(sparse_Portfolio_Cplex.w[i]^2/thePortfolioData.γ[sparse_Portfolio_Cplex.indices[i]] for i=1:size(sparse_Portfolio_Cplex.indices,1)))

        open("results/Table9ResultsOutput_withminreturn.csv","a") do fp
          println(fp, "Instance, cardinality, regularizer, method, solvetime, boundgap, boundgapSOCP, nodecount, number of iterations")
          println(fp, "rank"*string(rank), ",", thePortfolioData.k, ",",
                            "1/sqrt(n)", ",", "CPLEX MISOCP", ",",  sparse_Portfolio_Cplex.solveTime, ",", sparse_Portfolio_Cplex.boundGap, ",", sparse_Portfolio_Cplex.boundGapSOCP, ",", sparse_Portfolio_Cplex.nodeCount, ",", sparse_Portfolio_Cplex.numIters)
        end
    end
end


text:
struct FOne1 <: Function end; @inline (::FOne1)(x)=1.; f_one1 = FOne1()
struct FZero1 <: Function end; @inline (::FZero1)(x)=0.;  f_zero1 = FZero1()
struct FMOne1 <: Function end; @inline (::FMOne1)(x)=-1.;  f_mone1 = FMOne1()
struct FOne2 <: Function end; @inline (::FOne2)(x1,x2)=1.; f_one2 = FOne2()
struct FZero2 <: Function end; @inline (::FZero2)(x1,x2)=0.;  f_zero2 = FZero2()
struct FMOne2 <: Function end; @inline (::FMOne2)(x1,x2)=-1.;  f_mone2 = FMOne2()
struct FX1 <: Function end; @inline (::FX1)(x1,x2)=x1;  f_x1 = FX1()
struct FX2 <: Function end; @inline (::FX2)(x1,x2)=x2;  f_x2 = FX2()

const f_nargs_1 = [
    (
        :+,
        f_one1,
        f_zero1
    ),
    (
        :-,
        f_mone1,
        f_zero1
    ),
    (
        :inv,
        x -> -abs2(inv(x)),
        x -> -(2*inv(x))*(-abs2(inv(x)))
    ),
    (
        :abs2,
        x -> 2x,
        x -> 2
    ),
    (
        :exp,
        exp,
        exp
    ),
    (
        :exp2,
        x -> exp2(x)*0.69314718055994528622676398299518041312694549560546875,
        x -> exp2(x)*0.69314718055994528622676398299518041312694549560546875*0.69314718055994528622676398299518041312694549560546875
    ),
    (
        :exp10,
        x -> exp10(x)*2.30258509299404590109361379290930926799774169921875,
        x -> exp10(x)*2.30258509299404590109361379290930926799774169921875*2.30258509299404590109361379290930926799774169921875
    ),
    (
        :log,
        inv,
        x -> -abs2(inv(x))
    ),
    (
        :log2,
        x -> inv(x)/0.69314718055994528622676398299518041312694549560546875,
        x -> (-abs2(inv(x)))/0.69314718055994528622676398299518041312694549560546875
    ),
    (
        :log10,
        x -> inv(x)/2.30258509299404590109361379290930926799774169921875,
        x -> (-abs2(inv(x)))/2.30258509299404590109361379290930926799774169921875
    ),
    (
        :sin,
        cos,
        x -> -sin(x)
    ),
    (
        :cos,
        x -> -sin(x),
        x -> -cos(x)
    ),
    (
        :tan,
        x -> 1 + tan(x)^2,
        x -> 2*tan(x)*(1 + tan(x)^2)
    ),
    (
        :csc,
        x -> (-csc(x))*cot(x),
        x -> (-(-csc(x))*cot(x))*cot(x) + (-(1 + cot(x)^2))*(-csc(x))
    ),
    (
        :sec,
        x -> sec(x)*tan(x),
        x -> sec(x)*tan(x)*tan(x) + (1 + tan(x)^2)*sec(x)
    ),
    (
        :cot,
        x -> -(1 + cot(x)^2),
        x -> -2*cot(x)*(-(1 + cot(x)^2))
    ),
    (
        :atan,
        x -> inv(1 + x^2),
        x -> (-abs2(inv(1 + x^2)))*2x
    ),
    (
        :acot,
        x -> -inv(1 + x^2),
        x -> -(-abs2(inv(1 + x^2)))*2x
    ),
    (
        :sind,
        x -> 0.0174532925199432954743716805978692718781530857086181640625*cosd(x),
        x -> 0.0174532925199432954743716805978692718781530857086181640625*-0.0174532925199432954743716805978692718781530857086181640625*sind(x)
    ),
    (
        :cosd,
        x -> -0.0174532925199432954743716805978692718781530857086181640625*sind(x),
        x -> -0.0174532925199432954743716805978692718781530857086181640625*0.0174532925199432954743716805978692718781530857086181640625*cosd(x)
    ),
    (
        :tand,
        x -> 0.0174532925199432954743716805978692718781530857086181640625*(1 + tand(x)^2),
        x -> 0.0174532925199432954743716805978692718781530857086181640625*2*tand(x)*0.0174532925199432954743716805978692718781530857086181640625*(1 + tand(x)^2)
    ),
    (
        :cscd,
        x -> -0.0174532925199432954743716805978692718781530857086181640625*cscd(x)*cotd(x),
        x -> -0.0174532925199432954743716805978692718781530857086181640625*-0.0174532925199432954743716805978692718781530857086181640625*cscd(x)*cotd(x)*cotd(x) + -0.0174532925199432954743716805978692718781530857086181640625*(1 + cotd(x)^2)*-0.0174532925199432954743716805978692718781530857086181640625*cscd(x)
    ),
    (
        :secd,
        x -> 0.0174532925199432954743716805978692718781530857086181640625*secd(x)*tand(x),
        x -> 0.0174532925199432954743716805978692718781530857086181640625*0.0174532925199432954743716805978692718781530857086181640625*secd(x)*tand(x)*tand(x) + 0.0174532925199432954743716805978692718781530857086181640625*(1 + tand(x)^2)*0.0174532925199432954743716805978692718781530857086181640625*secd(x)
    ),
    (
        :cotd,
        x -> -0.0174532925199432954743716805978692718781530857086181640625*(1 + cotd(x)^2),
        x -> -0.0174532925199432954743716805978692718781530857086181640625*2*cotd(x)*-0.0174532925199432954743716805978692718781530857086181640625*(1 + cotd(x)^2)
    ),
    (
        :atand,
        x -> 57.29577951308232286464772187173366546630859375/(1 + x^2),
        x -> -57.29577951308232286464772187173366546630859375*2*x/(1 + x^2)^2
    ),
    (
        :acotd,
        x -> -57.29577951308232286464772187173366546630859375/(1 + x^2),
        x -> 57.29577951308232286464772187173366546630859375*2*x/(1 + x^2)^2
    ),
    (
        :sinh,
        cosh,
        sinh
    ),
    (
        :cosh,
        sinh,
        cosh
    ),
    (
        :tanh,
        x -> 1 - tanh(x)^2,
        x -> -2*tanh(x)*(1 - tanh(x)^2)
    ),
    (
        :csch,
        x -> (-coth(x))*csch(x),
        x -> csch(x)^2*csch(x) + (-coth(x))*csch(x)*(-coth(x))
    ),
    (
        :sech,
        x -> (-tanh(x))*sech(x),
        x -> (-(1 - tanh(x)^2))*sech(x) + (-tanh(x))*sech(x)*(-tanh(x))
    ),
    (
        :coth,
        x -> -csch(x)^2,
        x -> -2*csch(x)*(-coth(x))*csch(x)
    ),
    (
        :atanh,
        x -> inv(1 - x^2),
        x -> (-abs2(inv(1 - x^2)))*(-2x)
    ),
    (
        :acoth,
        x -> inv(1 - x^2),
        x -> (-abs2(inv(1 - x^2)))*(-2x)
    ),
    (
        :erfi,
        x -> 1.1283791670955125585606992899556644260883331298828125*exp(x^2),
        x -> 1.1283791670955125585606992899556644260883331298828125*exp(x^2)*2x
    ),
    (
        :loggamma,
        digamma,
        trigamma
    ),
    (
        :erfcinv,
        x -> -0.8862269254527579409597137782839126884937286376953125*exp(erfcinv(x)^2),
        x -> -0.8862269254527579409597137782839126884937286376953125*exp(erfcinv(x)^2)*2*erfcinv(x)*-0.8862269254527579409597137782839126884937286376953125*exp(erfcinv(x)^2)
    ),
    (
        :erfcx,
        x -> 2*x*erfcx(x) - 1.1283791670955125585606992899556644260883331298828125,
        x -> 2*erfcx(x) + (2*x*erfcx(x) - 1.1283791670955125585606992899556644260883331298828125)*2*x
    ),
    (
        :invdigamma,
        x -> inv(trigamma(invdigamma(x))),
        x -> (-abs2(inv(trigamma(invdigamma(x)))))*polygamma(2 , invdigamma(x))*inv(trigamma(invdigamma(x)))
    ),
    (
        :bessely1,
        x -> (bessely0(x) - bessely(2 , x))/2,
        x -> (-bessely1(x) + -(bessely(1 , x) - bessely(3 , x))/2)/2
    ),
    (
        :besselj1,
        x -> (besselj0(x) - besselj(2 , x))/2,
        x -> (-besselj1(x) + -(besselj(1 , x) - besselj(3 , x))/2)/2
    ),
    (
        :dawson,
        x -> 1 - 2*x*dawson(x),
        x -> -(2*dawson(x) + (1 - 2*x*dawson(x))*2*x)
    ),
    (
        :airyaiprime,
        x -> x*airyai(x),
        x -> airyai(x) + airyaiprime(x)*x
    ),
    (
        :erf,
        x -> 1.1283791670955125585606992899556644260883331298828125*exp(-x*x),
        x -> 1.1283791670955125585606992899556644260883331298828125*exp(-x*x)*(-2x)
    ),
    (
        :digamma,
        trigamma,
        x -> polygamma(2 , x)
    ),
    (
        :gamma,
        x -> digamma(x)*gamma(x),
        x -> trigamma(x)*gamma(x) + digamma(x)*gamma(x)*digamma(x)
    ),
    (
        :airyai,
        airyaiprime,
        x -> x*airyai(x)
    ),
    (
        :airybi,
        airybiprime,
        x -> x*airybi(x)
    ),
    (
        :erfinv,
        x -> 0.8862269254527579409597137782839126884937286376953125*exp(erfinv(x)^2),
        x -> 0.8862269254527579409597137782839126884937286376953125*exp(erfinv(x)^2)*2*erfinv(x)*0.8862269254527579409597137782839126884937286376953125*exp(erfinv(x)^2)
    ),
    (
        :bessely0,
        x -> -bessely1(x),
        x -> -(bessely0(x) - bessely(2 , x))/2
    ),
    (
        :erfc,
        x -> -1.1283791670955125585606992899556644260883331298828125*exp(-x*x),
        x -> -1.1283791670955125585606992899556644260883331298828125*exp(-x*x)*(-2x)
    ),
    (
        :trigamma,
        x -> polygamma(2 , x),
        x -> polygamma(3 , x)
    ),
    (
        :airybiprime,
        x -> x*airybi(x),
        x -> airybi(x) + airybiprime(x)*x
    ),
    (
        :besselj0,
        x -> -besselj1(x),
        x -> -(besselj0(x) - besselj(2 , x))/2
    )
]

const f_nargs_2 = [
    (
        :+,
        f_one2,
        f_one2,
        f_zero2,
        f_zero2,
        f_zero2
    ),
    (
        :-,
        f_one2,
        f_mone2,
        f_zero2,
        f_zero2,
        f_zero2
    ),
    (
        :*,
        f_x2,
        f_x1,
        f_zero2,
        f_one2,
        f_zero2
    ),
    (
        :^,
        ((x1,x2)  -> x2*x1^(x2 - 1)),
        ((x1,x2)  -> x1^x2*log(x1)),
        ((x1,x2)  -> x2*(x2 - 1)*x1^(x2 - 2)),
        ((x1,x2)  -> x2*x1^(x2  -  1)*log(x1) + x1^(x2 - 1)),
        ((x1,x2)  -> x1^x2*log(x1)*log(x1))),
    (
        :/,
        ((x1,x2)  -> 1/x2),
        ((x1,x2)  -> -x1/x2^2),
        f_zero2,
        ((x1,x2)  -> -1/x2^2),
        ((x1,x2)  -> 2x1/x2^3)
    ),
    (
        :beta,
        (x1,x2) -> beta(x1 , x2)*(digamma(x1) - digamma(x1 + x2)),
        (x1,x2) -> beta(x1 , x2)*(digamma(x2) - digamma(x1 + x2)),
        (x1,x2) -> beta(x1 , x2)*(digamma(x1) - digamma(x1 + x2))*(digamma(x1) - digamma(x1 + x2)) + (trigamma(x1) + -trigamma(x1 + x2))*beta(x1 , x2),
        (x1,x2) -> beta(x1 , x2)*(digamma(x2) - digamma(x1 + x2))*(digamma(x1) - digamma(x1 + x2)) + (-trigamma(x1 + x2))*beta(x1 , x2),
        (x1,x2) -> beta(x1 , x2)*(digamma(x2) - digamma(x1 + x2))*(digamma(x2) - digamma(x1 + x2)) + (trigamma(x2) + -trigamma(x1 + x2))*beta(x1 , x2)
    ),
    (
        :logbeta,
        (x1,x2) -> digamma(x1) - digamma(x1 + x2),
        (x1,x2) -> digamma(x2) - digamma(x1 + x2),
        (x1,x2) -> trigamma(x1) + -trigamma(x1 + x2),
        (x1,x2) -> -trigamma(x1 + x2),
        (x1,x2) -> trigamma(x2) + -trigamma(x1 + x2)
    )
]


text:
# # Gaussian Process for Event likelihoods
#
# ## Preliminary steps
#
# ### Loading necessary packages
using Plots
using AugmentedGaussianProcesses
using Distributions

# ## Creating some random data
n_data = 200
X = (rand(n_data) .- 0.5) * 40
r = 5.0
Y = rand.(NegativeBinomial.(r, AGP.logistic.(sin.(X))))
scatter(X, Y)

# ## Run GP model with negative binomial likelihood to learn p
kernel = SqExponentialKernel() ∘ ScaleTransform(1.0)
m_negbinomial = VGP(
    X, Y, kernel, NegBinomialLikelihood(r), AnalyticVI(); optimiser=false, verbose=2
)
@time train!(m_negbinomial, 20)

# ## Running the same model but with a Poisson likelihood
kernel = SqExponentialKernel() ∘ ScaleTransform(1.0)
m_poisson = VGP(
    X, Y, kernel, PoissonLikelihood(r), AnalyticVI(); optimiser=false, verbose=2
)
@time train!(m_poisson, 20)

# Prediction and plot function on a grid
# Create a grid and compute prediction on it
function compute_grid(model, n_grid=50)
    mins = -20
    maxs = 20
    x_grid = range(mins, maxs; length=n_grid) # Create a grid
    y_grid, sig_y_grid = proba_y(model, reshape(x_grid, :, 1)) # Predict the mean and variance on the grid
    return y_grid, sig_y_grid, x_grid
end

# Plot the data as a scatter plot
function plot_data(X, Y)
    return Plots.scatter(X, Y; alpha=0.33, msw=0.0, lab="", size=(800, 500))
end

function plot_model(model, X, Y, title=nothing)
    n_grid = 100
    y_grid, sig_y_grid, x_grid = compute_grid(model, n_grid)
    p = plot_data(X, Y)
    Plots.plot!(
        p,
        x_grid,
        y_grid;
        ribbon=2 * sqrt.(sig_y_grid), # Plot 2 std deviations
        title=title,
        color="red",
        lab="",
        linewidth=3.0,
    )
    return p
end;

# ## Comparison between the two likelihoods
Plots.plot(
    plot_model.(
        [m_negbinomial, m_poisson], Ref(X), Ref(Y), ["Neg. Binomial", "Poisson"]
    )...;
    layout=(1, 2),
)


text:
using Test
using Plots
using BenchmarkTools
using LinearAlgebra
using Random
using Roots

push!(LOAD_PATH, abspath(@__DIR__,"..","src"))

using SmoothedSpectralAbscissa ; const SSA = SmoothedSpectralAbscissa

##

"""
    lyap_simple!(A::AbstractMatrix)
returns the naive implementation of the Lyapunov objective function and its
derivative. Can be used for testing, plotting, etc
"""
function lyap_simple!(A::AbstractMatrix)
    function  myf(s::Float64)::Float64
        P = lyap(A - UniformScaling(s),one(A))
        tr(P)
    end
    function d_myf(s::Float64)::Float64
        sc,Id = UniformScaling(s), one(A)
        P =  lyap(A - sc,Id)
        Q = lyap(A' - sc,Id)
        -2.0tr(P*Q)
    end
    myf,d_myf
end

"""
    ssa_test(A::Matrix{Float64})

This is the easiest implementation
for testing and benchmarking.
"""
function ssa_test(A::Matrix{Float64} ; ssa_eps=nothing)
    myf,d_myf = lyap_simple!(A)
    ssa_eps = something(ssa_eps,SSA.default_eps_ssa(A))
    myg(s)= 1.0 / myf(s) - ssa_eps
    d_myg(s) = - d_myf(s) / myf(s)^2
    sa = SSA.spectral_abscissa(A)
    ssa_start = sa + 0.5*ssa_eps
    find_zero( (myg,d_myg),ssa_start, Roots.Newton())
end

function ssa_simple(A::AbstractMatrix, with_gradient::Bool,PQ,
        ssa_eps::Union{Nothing,Float64}=nothing)
    _ssa_eps = something(ssa_eps, SSA.default_eps_ssa(A))
    SSA.PQ_init!(A,PQ)
    _sa = SSA.spectral_abscissa(PQ)
    _start = _sa + 0.5*ssa_eps
    g(s)=SSA._g(s,PQ,_ssa_eps,_sa)
    x = range(_start, 10 ; length=200)
    y = g.(x)
    (x,y)
end

##


text:
module TestHLLE
    using ManifoldLearning
    using LinearAlgebra: eigvals
    using Test

    @testset "HLLE" begin

	k = 12
	d = 2
	X, L = swiss_roll()
	Y = transform(HLLE, X, k=k, d=d)

	@test outdim(Y) == d
	@test size(projection(Y), 2) == size(X, 2)
	@test neighbors(Y) == k
	@test length(eigvals(Y)) == d
    end
end


text:
"""
Base class for problem families.

It is an abstraction for problem parameterization (e.g by the
number of the search space dimensions) that allows
to instantiate `OptimizationProblem` for the concrete parameters.
"""
abstract type ProblemFamily{FS<:FitnessScheme} end

"""
Family of `FunctionBasedProblem` optimization problems
parameterized by the number of search space dimensions.
"""
struct FunctionBasedProblemFamily{F,FS<:FitnessScheme,FO} <: ProblemFamily{FS}
    objfunc::Function                     # Objective function
    name::String
    fitness_scheme::FS
    reserved_ss::RangePerDimSearchSpace   # search space for the first reserved dimensions
    range_per_dim::ParamBounds            # Default range per dimension
    opt_value::FO                         # optional optimal value, or nothing

    function FunctionBasedProblemFamily(
            objfunc::Function, name::String,
            fitness_scheme::FS, range::ParamBounds, opt_value::FO = nothing,
            reserved_ss::RangePerDimSearchSpace = ZERO_SEARCH_SPACE
    ) where {FS<:FitnessScheme, FO}
        if FO <: Number
            fitness_type(fitness_scheme) == FO ||
                throw(ArgumentError("Fitness type ($(fitness_type(fitness_scheme))) and opt_value type ($(FO)) do not match"))
        end
        new{FO, FS, FO}(objfunc, name, fitness_scheme,
                        reserved_ss, range, opt_value)
    end
end

objfunc(family::FunctionBasedProblemFamily) = family.objfunc

"""
    instantiate(family::FunctionBasedProblemFamily, ndim::Int)

Construct search space for `FunctionBasedProblem` with the given number of dimensions.
"""
function instantiate_search_space(family::FunctionBasedProblemFamily, ndim::Int)
    ndim >= numdims(family.reserved_ss) ||
        throw(ArgumentError("Cannot create $ndim-problem: number of dimensions less than reserved dimensions"))
    vcat(family.reserved_ss, symmetric_search_space(ndim - numdims(family.reserved_ss), family.range_per_dim))
end

"""
    instantiate(family::FunctionBasedProblemFamily, ndim::Int)

Construct `FunctionBasedProblem` with the given number of dimensions.
"""
instantiate(family::FunctionBasedProblemFamily, ndim::Int) =
    FunctionBasedProblem(family.objfunc, family.name, family.fitness_scheme,
                         instantiate_search_space(family, ndim), family.opt_value)

function instantiate(prob::OptimizationProblem, ndim::Int)
    numdims(prob) == ndim ||
        throw("Problem dimensions ($(numdims(prob))) and the requested dimension ($ndim) do not match")
    return prob # as is
end

@deprecate fixed_dim_problem instantiate

MinimizationProblemFamily(f::Function, name::String, range::ParamBounds, fmin::Float64) =
    FunctionBasedProblemFamily(f, name, MinimizingFitnessScheme, range, fmin)

MinimizationProblemFamily(f::Function, name::String, range::ParamBounds) =
    FunctionBasedProblemFamily(f, name, MinimizingFitnessScheme, range)

minimization_problem(f::Function, name::String, range::ParamBounds, ndim::Int) =
    instantiate(MinimizationProblemFamily(f, name, range), ndim)

minimization_problem(f::Function, range::ParamBounds, ndim::Int) =
    instantiate(MinimizationProblemFamily(f, "<unknown>", range), ndim)

minimization_problem(f::Function, name::String, range::ParamBounds, ndim::Int, fmin::Float64) =
    instantiate(MinimizationProblemFamily(f, name, range, fmin), ndim)

"""
    problem_set(ps::Dict{Any, FunctionBasedProblemFamily}, dims::Union{Int,Vector{Int}})

Construct a fixed-dimensional version of each problem from `ps`
for each dimension given in `dims`.

Returns a dictionary of problems.
"""
function problem_set(ps::Dict{Any, FunctionBasedProblemFamily}, dims::Vector{Int})
    next_free_index = 1
    # FIXME would using vector be more straightforward?
    result = Dict{Int, FunctionBasedProblem}()
    for d in dims
        for p in values(ps)
            result[next_free_index] = instantiate(p, d)
            next_free_index += 1
        end
    end
    return result
end

problem_set(ps::Dict{Any, FunctionBasedProblemFamily}, dim::Int) = problem_set(ps, [dim])


text:
push!(LOAD_PATH,"../src/")

using PetscWrap
using Documenter
using Literate

# Generate examples
example_src = joinpath(@__DIR__,"..","example")
example_dir = joinpath(@__DIR__,"src","example")
Sys.rm(example_dir; recursive=true, force=true)
Literate.markdown(joinpath(example_src, "linear_system.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
Literate.markdown(joinpath(example_src, "linear_system_fancy.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells

makedocs(;
    modules=[PetscWrap],
    authors="bmxam",
    sitename="PetscWrap.jl",
    clean=true,doctest=false,
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://github.com/bmxam/PetscWrap.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => Any[
            "example/linear_system.md",
        ],
        "Fancy examples" => Any[
            "example/linear_system_fancy.md",
        ],
        "API Reference" => Any[
            "api/init.md",
            "api/viewer.md",
            "api/vec.md",
            "api/mat.md",
            "api/ksp.md",
        ],
        "API fancy" => "api/fancy/fancy.md",
    ],
)

deploydocs(;
    repo="github.com/bmxam/PetscWrap.jl.git",
    push_preview = true
)