text:
# just the basics, Cube, Tetrahedron, Pyramid, AxesWire, Carpet

const Cube = (
   [
    Point3D(-0.5,  0.5, -0.5),
    Point3D(0.5,   0.5, -0.5),
    Point3D(0.5,  -0.5, -0.5),
    Point3D(-0.5, -0.5, -0.5),
    Point3D(-0.5,  0.5,  0.5),
    Point3D(0.5,   0.5,  0.5),
    Point3D(0.5,  -0.5,  0.5),
    Point3D(-0.5, -0.5,  0.5)
    ],
   [[1, 2, 3, 4],
    [2, 6, 7, 3],
    [6, 5, 8, 7],
    [5, 1, 4, 8],
    [1, 5, 6, 2],
    [4, 3, 7, 8]])

const Tetrahedron = (
    [
    Point3D(1, 1, 1),
    Point3D(-1, -1, 1),
    Point3D(-1, 1, -1),
    Point3D(1, -1, -1)
    ],
   [[1, 2, 3],
    [1, 2, 4],
    [2, 3, 4],
    [4, 1, 3]])

const Pyramid = ([
    Point3D(-1, -1, 0),
    Point3D(1, -1, 0),
    Point3D(1, 1, 0),
    Point3D(-1, 1, 0),
    Point3D(0, 0, 1)],
   [[1, 2, 3, 4],
    [1, 5, 2],
    [2, 5, 3],
    [3, 5, 4],
    [4, 5, 1]])

const Teapot = (
    [
    Point3D(0.0, -0.0, 0.488037),
    Point3D(0.003906, -0.042188, 0.476326),
    Point3D(0.003906, 0.042188, 0.476326),
    Point3D(0.010742, -0.0, 0.575333),
    Point3D(0.0125, -0.056251, 0.450561),
    Point3D(0.0125, 0.056251, 0.450561),
    Point3D(0.019531, -0.0, 0.413654),
    Point3D(0.021094, -0.042188, 0.424797),
    Point3D(0.021094, 0.042188, 0.424797),
    Point3D(0.025, -0.0, 0.413086),
    Point3D(0.03875, -0.19625, 0.488037),
    Point3D(0.03875, 0.19625, 0.488037),
    Point3D(0.039063, -0.0, 0.66803),
    Point3D(0.04866, -0.192034, 0.575333),
    Point3D(0.04866, 0.192034, 0.575333),
    Point3D(0.056768, -0.188584, 0.413654),
    Point3D(0.056768, 0.188584, 0.413654),
    Point3D(0.0625, -0.0, 0.358795),
    Point3D(0.074785, -0.180918, 0.66803),
    Point3D(0.074785, 0.180918, 0.66803),
    Point3D(0.079102, -0.0, 0.764481),
    Point3D(0.096406, -0.171719, 0.358795),
    Point3D(0.096406, 0.171719, 0.358795),
    Point3D(0.1, -0.0, 0.769043),
    Point3D(0.103906, -0.042188, 0.777779),
    Point3D(0.103906, 0.042188, 0.777779),
    Point3D(0.105469, -0.0, 0.32156),
    Point3D(0.111721, -0.165203, 0.764481),
    Point3D(0.111721, 0.165203, 0.764481),
    Point3D(0.1125, -0.056251, 0.796997),
    Point3D(0.1125, 0.056251, 0.796997),
    Point3D(0.121094, -0.042188, 0.816215),
    Point3D(0.121094, 0.042188, 0.816215),
    Point3D(0.125, -0.0, 0.300049),
    Point3D(0.125, -0.0, 0.824951),
    Point3D(0.125, -0.0, 0.863037),
    Point3D(0.136045, -0.154853, 0.32156),
    Point3D(0.136045, 0.154853, 0.32156),
    Point3D(0.137695, -0.0, 0.881027),
    Point3D(0.145, -0.355, 0.488037),
    Point3D(0.145, 0.355, 0.488037),
    Point3D(0.149219, -0.0, 0.887024),
    Point3D(0.15, -0.0, 0.863037),
    Point3D(0.152627, -0.347373, 0.575333),
    Point3D(0.152627, 0.347373, 0.575333),
    Point3D(0.154062, -0.147187, 0.300049),
    Point3D(0.154062, -0.147187, 0.863037),
    Point3D(0.154062, 0.147187, 0.300049),
    Point3D(0.154062, 0.147187, 0.863037),
    Point3D(0.154883, -0.0, 0.881027),
    Point3D(0.158867, -0.341133, 0.413654),
    Point3D(0.158867, 0.341133, 0.413654),
    Point3D(0.165774, -0.142205, 0.881027),
    Point3D(0.165774, 0.142204, 0.881027),
    Point3D(0.172734, -0.327266, 0.66803),
    Point3D(0.172734, 0.327266, 0.66803),
    Point3D(0.175, -0.0, 0.863037),
    Point3D(0.176404, -0.137682, 0.887024),
    Point3D(0.176404, 0.137681, 0.887024),
    Point3D(0.177125, -0.137375, 0.863037),
    Point3D(0.177125, 0.137375, 0.863037),
    Point3D(0.181629, -0.135459, 0.881027),
    Point3D(0.181629, 0.135458, 0.881027),
    Point3D(0.189375, -0.310625, 0.358795),
    Point3D(0.189375, 0.310625, 0.358795),
    Point3D(0.200187, -0.127562, 0.863037),
    Point3D(0.200187, 0.127562, 0.863037),
    Point3D(0.201162, -0.298838, 0.764481),
    Point3D(0.201162, 0.298838, 0.764481),
    Point3D(0.210938, -0.0, 0.884945),
    Point3D(0.219883, -0.280117, 0.32156),
    Point3D(0.219883, 0.280117, 0.32156),
    Point3D(0.23334, -0.113457, 0.884945),
    Point3D(0.23334, 0.113457, 0.884945),
    Point3D(0.23375, -0.26625, 0.300049),
    Point3D(0.23375, -0.26625, 0.863037),
    Point3D(0.23375, 0.26625, 0.300049),
    Point3D(0.23375, 0.26625, 0.863037),
    Point3D(0.242764, -0.257238, 0.881027),
    Point3D(0.242764, 0.257235, 0.881027),
    Point3D(0.250945, -0.249056, 0.887024),
    Point3D(0.250945, 0.249054, 0.887024),
    Point3D(0.2515, -0.2485, 0.863037),
    Point3D(0.2515, 0.2485, 0.863037),
    Point3D(0.254967, -0.245034, 0.881027),
    Point3D(0.254967, 0.245033, 0.881027),
    Point3D(0.26925, -0.23075, 0.863037),
    Point3D(0.26925, 0.23075, 0.863037),
    Point3D(0.29375, -0.0, 0.900055),
    Point3D(0.294766, -0.205234, 0.884945),
    Point3D(0.294766, 0.205234, 0.884945),
    Point3D(0.30375, -0.46125, 0.488037),
    Point3D(0.30375, 0.46125, 0.488037),
    Point3D(0.307966, -0.45134, 0.575333),
    Point3D(0.307966, 0.45134, 0.575333),
    Point3D(0.309734, -0.080953, 0.900055),
    Point3D(0.309734, 0.080953, 0.900055),
    Point3D(0.311416, -0.443233, 0.413654),
    Point3D(0.311416, 0.443233, 0.413654),
    Point3D(0.319082, -0.425215, 0.66803),
    Point3D(0.319082, 0.425215, 0.66803),
    Point3D(0.328281, -0.403594, 0.358795),
    Point3D(0.328281, 0.403594, 0.358795),
    Point3D(0.334797, -0.388279, 0.764481),
    Point3D(0.334797, 0.388279, 0.764481),
    Point3D(0.345146, -0.363955, 0.32156),
    Point3D(0.345146, 0.363955, 0.32156),
    Point3D(0.352812, -0.345937, 0.300049),
    Point3D(0.352812, -0.345937, 0.863037),
    Point3D(0.352812, 0.345937, 0.300049),
    Point3D(0.352812, 0.345937, 0.863037),
    Point3D(0.353562, -0.146438, 0.900055),
    Point3D(0.353562, 0.146438, 0.900055),
    Point3D(0.357795, -0.334229, 0.881027),
    Point3D(0.357795, 0.334224, 0.881027),
    Point3D(0.362318, -0.323599, 0.887024),
    Point3D(0.362318, 0.323595, 0.887024),
    Point3D(0.362625, -0.322875, 0.863037),
    Point3D(0.362625, 0.322875, 0.863037),
    Point3D(0.364542, -0.318373, 0.881027),
    Point3D(0.364542, 0.318371, 0.881027),
    Point3D(0.372437, -0.299813, 0.863037),
    Point3D(0.372437, 0.299813, 0.863037),
    Point3D(0.385938, -0.0, 0.915394),
    Point3D(0.386543, -0.266661, 0.884945),
    Point3D(0.386543, 0.266661, 0.884945),
    Point3D(0.394777, -0.044769, 0.915394),
    Point3D(0.394777, 0.044769, 0.915394),
    Point3D(0.414844, -0.0, 1.03776),
    Point3D(0.41875, -0.0, 1.00793),
    Point3D(0.419016, -0.080984, 0.915394),
    Point3D(0.419016, 0.080984, 0.915394),
    Point3D(0.419047, -0.190266, 0.900055),
    Point3D(0.419047, 0.190266, 0.900055),
    Point3D(0.421414, -0.033513, 1.03776),
    Point3D(0.421414, 0.033513, 1.03776),
    Point3D(0.425021, -0.03197, 1.00793),
    Point3D(0.425021, 0.03197, 1.00793),
    Point3D(0.43946, -0.06054, 1.03776),
    Point3D(0.43946, 0.06054, 1.03776),
    Point3D(0.442242, -0.057758, 1.00793),
    Point3D(0.442242, 0.057758, 1.00793),
    Point3D(0.45, -0.0, 0.937988),
    Point3D(0.450781, -0.0, 0.971153),
    Point3D(0.453875, -0.019625, 0.937988),
    Point3D(0.453875, 0.019625, 0.937988),
    Point3D(0.454586, -0.019348, 0.971153),
    Point3D(0.454586, 0.019348, 0.971153),
    Point3D(0.45523, -0.105222, 0.915394),
    Point3D(0.45523, 0.105222, 0.915394),
    Point3D(0.4645, -0.0355, 0.937988),
    Point3D(0.4645, 0.0355, 0.937988),
    Point3D(0.465028, -0.034972, 0.971153),
    Point3D(0.465028, 0.034972, 0.971153),
    Point3D(0.466487, -0.078587, 1.03776),
    Point3D(0.466487, 0.078587, 1.03776),
    Point3D(0.46803, -0.07498, 1.00793),
    Point3D(0.46803, 0.07498, 1.00793),
    Point3D(0.480375, -0.046124, 0.937988),
    Point3D(0.480375, 0.046124, 0.937988),
    Point3D(0.480652, -0.045414, 0.971153),
    Point3D(0.480652, 0.045414, 0.971153),
    Point3D(0.5, -0.0, 1.05005),
    Point3D(0.5, -0.049218, 0.971153),
    Point3D(0.5, -0.049999, 0.937988),
    Point3D(0.5, -0.08125, 1.00793),
    Point3D(0.5, -0.085157, 1.03776),
    Point3D(0.5, -0.114062, 0.915394),
    Point3D(0.5, -0.20625, 0.900055),
    Point3D(0.5, -0.289063, 0.884945),
    Point3D(0.5, -0.325001, 0.863037),
    Point3D(0.5, -0.34512, 0.881027),
    Point3D(0.5, -0.35, 0.863037),
    Point3D(0.5, -0.350785, 0.887024),
    Point3D(0.5, -0.362308, 0.881027),
    Point3D(0.5, -0.375, 0.300049),
    Point3D(0.5, -0.375, 0.863037),
    Point3D(0.5, -0.394531, 0.32156),
    Point3D(0.5, -0.420898, 0.764481),
    Point3D(0.5, -0.4375, 0.358795),
    Point3D(0.5, -0.460938, 0.66803),
    Point3D(0.5, -0.480469, 0.413654),
    Point3D(0.5, -0.489258, 0.575333),
    Point3D(0.5, -0.5, 0.488037),
    Point3D(0.5, 0.049218, 0.971153),
    Point3D(0.5, 0.049999, 0.937988),
    Point3D(0.5, 0.08125, 1.00793),
    Point3D(0.5, 0.085157, 1.03776),
    Point3D(0.5, 0.114062, 0.915394),
    Point3D(0.5, 0.20625, 0.900055),
    Point3D(0.5, 0.289063, 0.884945),
    Point3D(0.5, 0.325001, 0.863037),
    Point3D(0.5, 0.345118, 0.881027),
    Point3D(0.5, 0.35, 0.863037),
    Point3D(0.5, 0.35078, 0.887024),
    Point3D(0.5, 0.362303, 0.881027),
    Point3D(0.5, 0.375, 0.300049),
    Point3D(0.5, 0.375, 0.863037),
    Point3D(0.5, 0.394531, 0.32156),
    Point3D(0.5, 0.420898, 0.764481),
    Point3D(0.5, 0.4375, 0.358795),
    Point3D(0.5, 0.460938, 0.66803),
    Point3D(0.5, 0.480469, 0.413654),
    Point3D(0.5, 0.489258, 0.575333),
    Point3D(0.5, 0.5, 0.488037),
    Point3D(0.519348, -0.045414, 0.971153),
    Point3D(0.519348, 0.045414, 0.971153),
    Point3D(0.519625, -0.046124, 0.937988),
    Point3D(0.519625, 0.046124, 0.937988),
    Point3D(0.53197, -0.07498, 1.00793),
    Point3D(0.53197, 0.07498, 1.00793),
    Point3D(0.533513, -0.078587, 1.03776),
    Point3D(0.533513, 0.078587, 1.03776),
    Point3D(0.534972, -0.034972, 0.971153),
    Point3D(0.534972, 0.034972, 0.971153),
    Point3D(0.5355, -0.0355, 0.937988),
    Point3D(0.5355, 0.0355, 0.937988),
    Point3D(0.54477, -0.105222, 0.915394),
    Point3D(0.54477, 0.105222, 0.915394),
    Point3D(0.545414, -0.019348, 0.971153),
    Point3D(0.545414, 0.019348, 0.971153),
    Point3D(0.546125, -0.019625, 0.937988),
    Point3D(0.546125, 0.019625, 0.937988),
    Point3D(0.549219, -0.0, 0.971153),
    Point3D(0.55, -0.0, 0.937988),
    Point3D(0.557758, -0.057758, 1.00793),
    Point3D(0.557758, 0.057758, 1.00793),
    Point3D(0.56054, -0.06054, 1.03776),
    Point3D(0.56054, 0.06054, 1.03776),
    Point3D(0.574979, -0.03197, 1.00793),
    Point3D(0.574979, 0.03197, 1.00793),
    Point3D(0.578586, -0.033513, 1.03776),
    Point3D(0.578586, 0.033513, 1.03776),
    Point3D(0.580953, -0.190266, 0.900055),
    Point3D(0.580953, 0.190266, 0.900055),
    Point3D(0.580984, -0.080984, 0.915394),
    Point3D(0.580984, 0.080984, 0.915394),
    Point3D(0.58125, -0.0, 1.00793),
    Point3D(0.585156, -0.0, 1.03776),
    Point3D(0.605223, -0.044769, 0.915394),
    Point3D(0.605223, 0.044769, 0.915394),
    Point3D(0.613457, -0.266661, 0.884945),
    Point3D(0.613457, 0.266661, 0.884945),
    Point3D(0.614062, -0.0, 0.915394),
    Point3D(0.627562, -0.299813, 0.863037),
    Point3D(0.627562, 0.299813, 0.863037),
    Point3D(0.635459, -0.318373, 0.881027),
    Point3D(0.635459, 0.318371, 0.881027),
    Point3D(0.637375, -0.322875, 0.863037),
    Point3D(0.637375, 0.322875, 0.863037),
    Point3D(0.637682, -0.323599, 0.887024),
    Point3D(0.637682, 0.323595, 0.887024),
    Point3D(0.642205, -0.334229, 0.881027),
    Point3D(0.642205, 0.334224, 0.881027),
    Point3D(0.646437, -0.146438, 0.900055),
    Point3D(0.646437, 0.146438, 0.900055),
    Point3D(0.647187, -0.345937, 0.300049),
    Point3D(0.647187, -0.345937, 0.863037),
    Point3D(0.647187, 0.345937, 0.300049),
    Point3D(0.647187, 0.345937, 0.863037),
    Point3D(0.654853, -0.363955, 0.32156),
    Point3D(0.654853, 0.363955, 0.32156),
    Point3D(0.665203, -0.388279, 0.764481),
    Point3D(0.665203, 0.388279, 0.764481),
    Point3D(0.671719, -0.403594, 0.358795),
    Point3D(0.671719, 0.403594, 0.358795),
    Point3D(0.680918, -0.425215, 0.66803),
    Point3D(0.680918, 0.425215, 0.66803),
    Point3D(0.688584, -0.443233, 0.413654),
    Point3D(0.688584, 0.443233, 0.413654),
    Point3D(0.690266, -0.080953, 0.900055),
    Point3D(0.690266, 0.080953, 0.900055),
    Point3D(0.692034, -0.45134, 0.575333),
    Point3D(0.692034, 0.45134, 0.575333),
    Point3D(0.69625, -0.46125, 0.488037),
    Point3D(0.69625, 0.46125, 0.488037),
    Point3D(0.705234, -0.205234, 0.884945),
    Point3D(0.705234, 0.205234, 0.884945),
    Point3D(0.70625, -0.0, 0.900055),
    Point3D(0.73075, -0.23075, 0.863037),
    Point3D(0.73075, 0.23075, 0.863037),
    Point3D(0.745033, -0.245034, 0.881027),
    Point3D(0.745033, 0.245033, 0.881027),
    Point3D(0.7485, -0.2485, 0.863037),
    Point3D(0.7485, 0.2485, 0.863037),
    Point3D(0.749055, -0.249056, 0.887024),
    Point3D(0.749055, 0.249054, 0.887024),
    Point3D(0.757236, -0.257238, 0.881027),
    Point3D(0.757236, 0.257235, 0.881027),
    Point3D(0.76625, -0.26625, 0.300049),
    Point3D(0.76625, -0.26625, 0.863037),
    Point3D(0.76625, 0.26625, 0.300049),
    Point3D(0.76625, 0.26625, 0.863037),
    Point3D(0.76666, -0.113457, 0.884945),
    Point3D(0.76666, 0.113457, 0.884945),
    Point3D(0.780117, -0.280117, 0.32156),
    Point3D(0.780117, 0.280117, 0.32156),
    Point3D(0.789062, -0.0, 0.884945),
    Point3D(0.798838, -0.298838, 0.764481),
    Point3D(0.798838, 0.298838, 0.764481),
    Point3D(0.799812, -0.127562, 0.863037),
    Point3D(0.799812, 0.127562, 0.863037),
    Point3D(0.810625, -0.310625, 0.358795),
    Point3D(0.810625, 0.310625, 0.358795),
    Point3D(0.818371, -0.135459, 0.881027),
    Point3D(0.818371, 0.135458, 0.881027),
    Point3D(0.822875, -0.137375, 0.863037),
    Point3D(0.822875, 0.137375, 0.863037),
    Point3D(0.823596, -0.137682, 0.887024),
    Point3D(0.823596, 0.137681, 0.887024),
    Point3D(0.825, -0.0, 0.863037),
    Point3D(0.827266, -0.327266, 0.66803),
    Point3D(0.827266, 0.327266, 0.66803),
    Point3D(0.834226, -0.142205, 0.881027),
    Point3D(0.834226, 0.142204, 0.881027),
    Point3D(0.841133, -0.341133, 0.413654),
    Point3D(0.841133, 0.341133, 0.413654),
    Point3D(0.845117, -0.0, 0.881027),
    Point3D(0.845937, -0.147187, 0.300049),
    Point3D(0.845937, -0.147187, 0.863037),
    Point3D(0.845937, 0.147187, 0.300049),
    Point3D(0.845937, 0.147187, 0.863037),
    Point3D(0.847373, -0.347373, 0.575333),
    Point3D(0.847373, 0.347373, 0.575333),
    Point3D(0.85, -0.0, 0.863037),
    Point3D(0.850781, -0.0, 0.887024),
    Point3D(0.855, -0.355, 0.488037),
    Point3D(0.855, 0.355, 0.488037),
    Point3D(0.862305, -0.0, 0.881027),
    Point3D(0.863955, -0.154853, 0.32156),
    Point3D(0.863955, 0.154853, 0.32156),
    Point3D(0.875, -0.0, 0.300049),
    Point3D(0.875, -0.0, 0.863037),
    Point3D(0.888279, -0.165203, 0.764481),
    Point3D(0.888279, 0.165203, 0.764481),
    Point3D(0.894531, -0.0, 0.32156),
    Point3D(0.903594, -0.171719, 0.358795),
    Point3D(0.903594, 0.171719, 0.358795),
    Point3D(0.920898, -0.0, 0.764481),
    Point3D(0.925, -0.0, 0.413086),
    Point3D(0.925, -0.0, 0.618896),
    Point3D(0.925, -0.092813, 0.445244),
    Point3D(0.925, -0.092813, 0.586738),
    Point3D(0.925, -0.123751, 0.515991),
    Point3D(0.925, 0.092813, 0.445244),
    Point3D(0.925, 0.092813, 0.586738),
    Point3D(0.925, 0.123751, 0.515991),
    Point3D(0.925215, -0.180918, 0.66803),
    Point3D(0.925215, 0.180918, 0.66803),
    Point3D(0.9375, -0.0, 0.358795),
    Point3D(0.943232, -0.188584, 0.413654),
    Point3D(0.943232, 0.188584, 0.413654),
    Point3D(0.95134, -0.192034, 0.575333),
    Point3D(0.95134, 0.192034, 0.575333),
    Point3D(0.960938, -0.0, 0.66803),
    Point3D(0.96125, -0.19625, 0.488037),
    Point3D(0.96125, 0.19625, 0.488037),
    Point3D(0.980469, -0.0, 0.413654),
    Point3D(0.989258, -0.0, 0.575333),
    Point3D(1.0, -0.0, 0.488037),
    Point3D(1.04492, -0.0, 0.646503),
    Point3D(1.05408, -0.083804, 0.622637),
    Point3D(1.05408, 0.083804, 0.622637),
    Point3D(1.07422, -0.111739, 0.570131),
    Point3D(1.07422, 0.111739, 0.570131),
    Point3D(1.09436, -0.083804, 0.517625),
    Point3D(1.09436, 0.083804, 0.517625),
    Point3D(1.09687, -0.0, 0.71286),
    Point3D(1.10352, -0.0, 0.493759),
    Point3D(1.10859, -0.063985, 0.698941),
    Point3D(1.10859, 0.063985, 0.698941),
    Point3D(1.12539, -0.0, 0.79327),
    Point3D(1.13437, -0.085313, 0.66832),
    Point3D(1.13437, 0.085313, 0.66832),
    Point3D(1.13967, -0.044165, 0.788218),
    Point3D(1.13967, 0.044165, 0.788218),
    Point3D(1.16016, -0.063985, 0.637698),
    Point3D(1.16016, 0.063985, 0.637698),
    Point3D(1.17109, -0.058887, 0.777105),
    Point3D(1.17109, 0.058887, 0.777105),
    Point3D(1.17188, -0.0, 0.623779),
    Point3D(1.175, -0.0, 0.863037),
    Point3D(1.19297, -0.0, 0.8732),
    Point3D(1.19844, -0.035156, 0.863037),
    Point3D(1.19844, 0.035156, 0.863037),
    Point3D(1.2, -0.0, 0.863037),
    Point3D(1.20251, -0.044165, 0.765992),
    Point3D(1.20251, 0.044165, 0.765992),
    Point3D(1.20625, -0.0, 0.876587),
    Point3D(1.21016, -0.0, 0.8732),
    Point3D(1.21563, -0.021094, 0.863037),
    Point3D(1.21563, 0.021094, 0.863037),
    Point3D(1.2168, -0.0, 0.76094),
    Point3D(1.21807, -0.032959, 0.873741),
    Point3D(1.21807, 0.032959, 0.873741),
    Point3D(1.22981, -0.028125, 0.87746),
    Point3D(1.22981, 0.028125, 0.87746),
    Point3D(1.23016, -0.023291, 0.873967),
    Point3D(1.23016, 0.023291, 0.873967),
    Point3D(1.25, -0.028125, 0.863037),
    Point3D(1.25, -0.046875, 0.863037),
    Point3D(1.25, 0.028125, 0.863037),
    Point3D(1.25, 0.046875, 0.863037),
    Point3D(1.27329, -0.043945, 0.874933),
    Point3D(1.27329, 0.043945, 0.874933),
    Point3D(1.27417, -0.031055, 0.875654),
    Point3D(1.27417, 0.031055, 0.875654),
    Point3D(1.28164, -0.0375, 0.879379),
    Point3D(1.28164, 0.0375, 0.879379),
    Point3D(1.28437, -0.021094, 0.863037),
    Point3D(1.28437, 0.021094, 0.863037),
    Point3D(1.3, -0.0, 0.863037),
    Point3D(1.30156, -0.035156, 0.863037),
    Point3D(1.30156, 0.035156, 0.863037),
    Point3D(1.31818, -0.023291, 0.877342),
    Point3D(1.31818, 0.023291, 0.877342),
    Point3D(1.325, -0.0, 0.863037),
    Point3D(1.32851, -0.032959, 0.876125),
    Point3D(1.32851, 0.032959, 0.876125),
    Point3D(1.33347, -0.028125, 0.881299),
    Point3D(1.33347, 0.028125, 0.881299),
    Point3D(1.33818, -0.0, 0.878109),
    Point3D(1.35361, -0.0, 0.876667),
    Point3D(1.35703, -0.0, 0.882171),
    Point3D(-0.016797, -0.0, 0.768166),
    Point3D(-0.020276, -0.042188, 0.776765),
    Point3D(-0.020276, 0.042188, 0.776765),
    Point3D(-0.02793, -0.056251, 0.795683),
    Point3D(-0.02793, 0.056251, 0.795683),
    Point3D(-0.035583, -0.042188, 0.814601),
    Point3D(-0.035583, 0.042188, 0.814601),
    Point3D(-0.039063, -0.0, 0.8232),
    Point3D(-0.080078, -0.0, 0.538757),
    Point3D(-0.082733, -0.042188, 0.529498),
    Point3D(-0.082733, 0.042188, 0.529498),
    Point3D(-0.088574, -0.056251, 0.509128),
    Point3D(-0.088574, 0.056251, 0.509128),
    Point3D(-0.094415, -0.042188, 0.488759),
    Point3D(-0.094415, 0.042188, 0.488759),
    Point3D(-0.09707, -0.0, 0.4795),
    Point3D(-0.103125, -0.0, 0.762024),
    Point3D(-0.111426, -0.042188, 0.769668),
    Point3D(-0.111426, 0.042188, 0.769668),
    Point3D(-0.129688, -0.056251, 0.786484),
    Point3D(-0.129688, 0.056251, 0.786484),
    Point3D(-0.134375, -0.0, 0.600098),
    Point3D(-0.141943, -0.042188, 0.592611),
    Point3D(-0.141943, 0.042188, 0.592611),
    Point3D(-0.147949, -0.042188, 0.8033),
    Point3D(-0.147949, 0.042188, 0.8033),
    Point3D(-0.15625, -0.0, 0.810943),
    Point3D(-0.156641, -0.0, 0.745354),
    Point3D(-0.158594, -0.056251, 0.576141),
    Point3D(-0.158594, 0.056251, 0.576141),
    Point3D(-0.165234, -0.0, 0.661621),
    Point3D(-0.167566, -0.042188, 0.750404),
    Point3D(-0.167566, 0.042188, 0.750404),
    Point3D(-0.175, -0.0, 0.712891),
    Point3D(-0.175244, -0.042188, 0.559671),
    Point3D(-0.175244, 0.042188, 0.559671),
    Point3D(-0.175885, -0.042188, 0.656723),
    Point3D(-0.175885, 0.042188, 0.656723),
    Point3D(-0.182812, -0.0, 0.552185),
    Point3D(-0.186719, -0.042188, 0.712891),
    Point3D(-0.186719, 0.042188, 0.712891),
    Point3D(-0.191602, -0.056251, 0.761515),
    Point3D(-0.191602, 0.056251, 0.761515),
    Point3D(-0.199316, -0.056251, 0.645947),
    Point3D(-0.199316, 0.056251, 0.645947),
    Point3D(-0.2125, -0.056251, 0.712891),
    Point3D(-0.2125, 0.056251, 0.712891),
    Point3D(-0.215637, -0.042188, 0.772625),
    Point3D(-0.215637, 0.042188, 0.772625),
    Point3D(-0.222748, -0.042188, 0.63517),
    Point3D(-0.222748, 0.042188, 0.63517),
    Point3D(-0.226562, -0.0, 0.777676),
    Point3D(-0.233398, -0.0, 0.630272),
    Point3D(-0.238281, -0.042188, 0.712891),
    Point3D(-0.238281, 0.042188, 0.712891),
    Point3D(-0.25, -0.0, 0.712891)],
    Array{Int64,1}[
        [325, 318, 305, 307],
        [307, 305, 282, 284],
        [284, 282, 247, 249],
        [249, 247, 172, 173],
        [318, 326, 309, 305],
        [305, 309, 286, 282],
        [282, 286, 251, 247],
        [247, 251, 174, 172],
        [326, 329, 314, 309],
        [309, 314, 288, 286],
        [286, 288, 253, 251],
        [251, 253, 175, 174],
        [329, 333, 320, 314],
        [314, 320, 291, 288],
        [288, 291, 258, 253],
        [253, 258, 177, 175],
        [173, 172, 120, 118],
        [118, 120, 85, 83],
        [83, 85, 62, 60],
        [60, 62, 50, 43],
        [172, 174, 116, 120],
        [120, 116, 81, 85],
        [85, 81, 58, 62],
        [62, 58, 42, 50],
        [174, 175, 114, 116],
        [116, 114, 79, 81],
        [81, 79, 53, 58],
        [58, 53, 39, 42],
        [175, 177, 109, 114],
        [114, 109, 76, 79],
        [79, 76, 47, 53],
        [53, 47, 36, 39],
        [43, 50, 63, 61],
        [61, 63, 86, 84],
        [84, 86, 121, 119],
        [119, 121, 193, 194],
        [50, 42, 59, 63],
        [63, 59, 82, 86],
        [86, 82, 117, 121],
        [121, 117, 195, 193],
        [42, 39, 54, 59],
        [59, 54, 80, 82],
        [82, 80, 115, 117],
        [117, 115, 196, 195],
        [39, 36, 49, 54],
        [54, 49, 78, 80],
        [80, 78, 111, 115],
        [115, 111, 198, 196],
        [194, 193, 248, 250],
        [250, 248, 283, 285],
        [285, 283, 306, 308],
        [308, 306, 318, 325],
        [193, 195, 252, 248],
        [248, 252, 287, 283],
        [283, 287, 310, 306],
        [306, 310, 326, 318],
        [195, 196, 254, 252],
        [252, 254, 289, 287],
        [287, 289, 315, 310],
        [310, 315, 329, 326],
        [196, 198, 260, 254],
        [254, 260, 293, 289],
        [289, 293, 322, 315],
        [315, 322, 333, 329],
        [333, 339, 334, 320],
        [320, 334, 299, 291],
        [291, 299, 263, 258],
        [258, 263, 179, 177],
        [339, 355, 348, 334],
        [334, 348, 312, 299],
        [299, 312, 267, 263],
        [263, 267, 181, 179],
        [355, 359, 353, 348],
        [348, 353, 323, 312],
        [312, 323, 273, 267],
        [267, 273, 183, 181],
        [359, 360, 356, 353],
        [353, 356, 327, 323],
        [323, 327, 275, 273],
        [273, 275, 184, 183],
        [177, 179, 104, 109],
        [109, 104, 68, 76],
        [76, 68, 28, 47],
        [47, 28, 21, 36],
        [179, 181, 100, 104],
        [104, 100, 55, 68],
        [68, 55, 19, 28],
        [28, 19, 13, 21],
        [181, 183, 94, 100],
        [100, 94, 44, 55],
        [55, 44, 14, 19],
        [19, 14, 4, 13],
        [183, 184, 92, 94],
        [94, 92, 40, 44],
        [44, 40, 11, 14],
        [14, 11, 1, 4],
        [36, 21, 29, 49],
        [49, 29, 69, 78],
        [78, 69, 105, 111],
        [111, 105, 200, 198],
        [21, 13, 20, 29],
        [29, 20, 56, 69],
        [69, 56, 101, 105],
        [105, 101, 202, 200],
        [13, 4, 15, 20],
        [20, 15, 45, 56],
        [56, 45, 95, 101],
        [101, 95, 204, 202],
        [4, 1, 12, 15],
        [15, 12, 41, 45],
        [45, 41, 93, 95],
        [95, 93, 205, 204],
        [198, 200, 264, 260],
        [260, 264, 300, 293],
        [293, 300, 335, 322],
        [322, 335, 339, 333],
        [200, 202, 268, 264],
        [264, 268, 313, 300],
        [300, 313, 349, 335],
        [335, 349, 355, 339],
        [202, 204, 274, 268],
        [268, 274, 324, 313],
        [313, 324, 354, 349],
        [349, 354, 359, 355],
        [204, 205, 276, 274],
        [274, 276, 328, 324],
        [324, 328, 357, 354],
        [354, 357, 360, 359],
        [360, 358, 351, 356],
        [356, 351, 316, 327],
        [327, 316, 269, 275],
        [275, 269, 182, 184],
        [358, 350, 337, 351],
        [351, 337, 303, 316],
        [316, 303, 265, 269],
        [269, 265, 180, 182],
        [350, 336, 330, 337],
        [337, 330, 296, 303],
        [303, 296, 261, 265],
        [265, 261, 178, 180],
        [336, 332, 319, 330],
        [330, 319, 290, 296],
        [296, 290, 257, 261],
        [261, 257, 176, 178],
        [184, 182, 98, 92],
        [92, 98, 51, 40],
        [40, 51, 16, 11],
        [11, 16, 7, 1],
        [182, 180, 102, 98],
        [98, 102, 64, 51],
        [51, 64, 22, 16],
        [16, 22, 18, 7],
        [180, 178, 106, 102],
        [102, 106, 71, 64],
        [64, 71, 37, 22],
        [22, 37, 27, 18],
        [178, 176, 108, 106],
        [106, 108, 75, 71],
        [71, 75, 46, 37],
        [37, 46, 34, 27],
        [1, 7, 17, 12],
        [12, 17, 52, 41],
        [41, 52, 99, 93],
        [93, 99, 203, 205],
        [7, 18, 23, 17],
        [17, 23, 65, 52],
        [52, 65, 103, 99],
        [99, 103, 201, 203],
        [18, 27, 38, 23],
        [23, 38, 72, 65],
        [65, 72, 107, 103],
        [103, 107, 199, 201],
        [27, 34, 48, 38],
        [38, 48, 77, 72],
        [72, 77, 110, 107],
        [107, 110, 197, 199],
        [205, 203, 270, 276],
        [276, 270, 317, 328],
        [328, 317, 352, 357],
        [357, 352, 358, 360],
        [203, 201, 266, 270],
        [270, 266, 304, 317],
        [317, 304, 338, 352],
        [352, 338, 350, 358],
        [201, 199, 262, 266],
        [266, 262, 297, 304],
        [304, 297, 331, 338],
        [338, 331, 336, 350],
        [199, 197, 259, 262],
        [262, 259, 292, 297],
        [297, 292, 321, 331],
        [331, 321, 332, 336],
        [24, 425, 426, 25],
        [25, 426, 428, 30],
        [30, 428, 430, 32],
        [32, 430, 432, 35],
        [425, 441, 442, 426],
        [426, 442, 444, 428],
        [428, 444, 449, 430],
        [430, 449, 451, 432],
        [441, 452, 456, 442],
        [442, 456, 466, 444],
        [444, 466, 472, 449],
        [449, 472, 476, 451],
        [452, 458, 464, 456],
        [456, 464, 470, 466],
        [466, 470, 478, 472],
        [472, 478, 480, 476],
        [35, 432, 431, 33],
        [33, 431, 429, 31],
        [31, 429, 427, 26],
        [26, 427, 425, 24],
        [432, 451, 450, 431],
        [431, 450, 445, 429],
        [429, 445, 443, 427],
        [427, 443, 441, 425],
        [451, 476, 473, 450],
        [450, 473, 467, 445],
        [445, 467, 457, 443],
        [443, 457, 452, 441],
        [476, 480, 479, 473],
        [473, 479, 471, 467],
        [467, 471, 465, 457],
        [457, 465, 458, 452],
        [458, 455, 461, 464],
        [464, 461, 468, 470],
        [470, 468, 474, 478],
        [478, 474, 477, 480],
        [455, 446, 447, 461],
        [461, 447, 453, 468],
        [468, 453, 459, 474],
        [474, 459, 463, 477],
        [446, 433, 434, 447],
        [447, 434, 436, 453],
        [453, 436, 438, 459],
        [459, 438, 440, 463],
        [433, 1, 2, 434],
        [434, 2, 5, 436],
        [436, 5, 8, 438],
        [438, 8, 10, 440],
        [480, 477, 475, 479],
        [479, 475, 469, 471],
        [471, 469, 462, 465],
        [465, 462, 455, 458],
        [477, 463, 460, 475],
        [475, 460, 454, 469],
        [469, 454, 448, 462],
        [462, 448, 446, 455],
        [463, 440, 439, 460],
        [460, 439, 437, 454],
        [454, 437, 435, 448],
        [448, 435, 433, 446],
        [440, 10, 9, 439],
        [439, 9, 6, 437],
        [437, 6, 3, 435],
        [435, 3, 1, 433],
        [341, 361, 362, 343],
        [343, 362, 364, 344],
        [344, 364, 366, 342],
        [342, 366, 369, 340],
        [361, 368, 370, 362],
        [362, 370, 373, 364],
        [364, 373, 377, 366],
        [366, 377, 381, 369],
        [368, 372, 375, 370],
        [370, 375, 379, 373],
        [373, 379, 387, 377],
        [377, 387, 393, 381],
        [372, 382, 384, 375],
        [375, 384, 401, 379],
        [379, 401, 413, 387],
        [387, 413, 417, 393],
        [340, 369, 367, 345],
        [345, 367, 365, 347],
        [347, 365, 363, 346],
        [346, 363, 361, 341],
        [369, 381, 378, 367],
        [367, 378, 374, 365],
        [365, 374, 371, 363],
        [363, 371, 368, 361],
        [381, 393, 388, 378],
        [378, 388, 380, 374],
        [374, 380, 376, 371],
        [371, 376, 372, 368],
        [393, 417, 414, 388],
        [388, 414, 403, 380],
        [380, 403, 385, 376],
        [376, 385, 382, 372],
        [382, 383, 394, 384],
        [384, 394, 404, 401],
        [401, 404, 418, 413],
        [413, 418, 423, 417],
        [383, 389, 396, 394],
        [394, 396, 408, 404],
        [404, 408, 420, 418],
        [418, 420, 424, 423],
        [389, 390, 398, 396],
        [396, 398, 406, 408],
        [408, 406, 415, 420],
        [420, 415, 422, 424],
        [390, 386, 391, 398],
        [398, 391, 400, 406],
        [406, 400, 410, 415],
        [415, 410, 412, 422],
        [417, 423, 419, 414],
        [414, 419, 405, 403],
        [403, 405, 395, 385],
        [385, 395, 383, 382],
        [423, 424, 421, 419],
        [419, 421, 409, 405],
        [405, 409, 397, 395],
        [395, 397, 389, 383],
        [424, 422, 416, 421],
        [421, 416, 407, 409],
        [409, 407, 399, 397],
        [397, 399, 390, 389],
        [422, 412, 411, 416],
        [416, 411, 402, 407],
        [407, 402, 392, 399],
        [399, 392, 386, 390],
        [163, 239, 232],
        [163, 232, 228],
        [163, 228, 212],
        [163, 212, 167],
        [239, 238, 230, 232],
        [232, 230, 226, 228],
        [228, 226, 210, 212],
        [212, 210, 166, 167],
        [238, 224, 220, 230],
        [230, 220, 214, 226],
        [226, 214, 206, 210],
        [210, 206, 164, 166],
        [224, 225, 222, 220],
        [220, 222, 216, 214],
        [214, 216, 208, 206],
        [206, 208, 165, 164],
        [163, 167, 155],
        [163, 155, 139],
        [163, 139, 135],
        [163, 135, 129],
        [167, 166, 157, 155],
        [155, 157, 141, 139],
        [139, 141, 137, 135],
        [135, 137, 130, 129],
        [166, 164, 161, 157],
        [157, 161, 153, 141],
        [141, 153, 147, 137],
        [137, 147, 144, 130],
        [164, 165, 159, 161],
        [161, 159, 151, 153],
        [153, 151, 145, 147],
        [147, 145, 143, 144],
        [163, 129, 136],
        [163, 136, 140],
        [163, 140, 156],
        [163, 156, 188],
        [129, 130, 138, 136],
        [136, 138, 142, 140],
        [140, 142, 158, 156],
        [156, 158, 187, 188],
        [130, 144, 148, 138],
        [138, 148, 154, 142],
        [142, 154, 162, 158],
        [158, 162, 185, 187],
        [144, 143, 146, 148],
        [148, 146, 152, 154],
        [154, 152, 160, 162],
        [162, 160, 186, 185],
        [163, 188, 213],
        [163, 213, 229],
        [163, 229, 233],
        [163, 233, 239],
        [188, 187, 211, 213],
        [213, 211, 227, 229],
        [229, 227, 231, 233],
        [233, 231, 238, 239],
        [187, 185, 207, 211],
        [211, 207, 215, 227],
        [227, 215, 221, 231],
        [231, 221, 224, 238], [185, 186, 209, 207], [207, 209, 217, 215], [215, 217, 223, 221], [221, 223, 225, 224], [225, 244, 240, 222], [222, 240, 236, 216], [216, 236, 218, 208], [208, 218, 168, 165], [244, 279, 271, 240], [240, 271, 255, 236], [236, 255, 234, 218], [218, 234, 169, 168], [279, 298, 294, 271], [271, 294, 277, 255], [255, 277, 242, 234], [234, 242, 170, 169], [298, 311, 301, 294], [294, 301, 280, 277], [277, 280, 245, 242], [242, 245, 171, 170], [165, 168, 149, 159], [159, 149, 131, 151], [151, 131, 127, 145], [145, 127, 124, 143], [168, 169, 133, 149], [149, 133, 112, 131], [131, 112, 96, 127], [127, 96, 89, 124], [169, 170, 125, 133], [133, 125, 90, 112], [112, 90, 73, 96], [96, 73, 70, 89], [170, 171, 122, 125], [125, 122, 87, 90], [90, 87, 66, 73], [73, 66, 57, 70], [143, 124, 128, 146], [146, 128, 132, 152], [152, 132, 150, 160], [160, 150, 189, 186], [124, 89, 97, 128], [128, 97, 113, 132], [132, 113, 134, 150], [150, 134, 190, 189], [89, 70, 74, 97],
    [97, 74, 91, 113],
    [113, 91, 126, 134], [134, 126, 191, 190], [70, 57, 67, 74], [74, 67, 88, 91], [91, 88, 123, 126], [126, 123, 192, 191], [186, 189, 219, 209], [209, 219, 237, 217], [217, 237, 241, 223], [223, 241, 244, 225], [189, 190, 235, 219], [219, 235, 256, 237], [237, 256, 272, 241], [241, 272, 279, 244], [190, 191, 243, 235], [235, 243, 278, 256], [256, 278, 295, 272], [272, 295, 298, 279], [191, 192, 246, 243], [243, 246, 281, 278], [278, 281, 302, 295], [295, 302, 311, 298]])


text:
# Raw memory management

export Mem

module Mem

using ..CUDAdrv
import ..@apicall, ..CuStream_t, ..CuDevice_t

## buffer type

struct Buffer
    ptr::CuPtr{Cvoid}
    bytesize::Int

    ctx::CuContext
end

Base.unsafe_convert(::Type{<:Ptr}, buf::Buffer) =
    throw(ArgumentError("cannot take the CPU address of a GPU buffer"))
Base.unsafe_convert(::Type{CuPtr{T}}, buf::Buffer) where {T} = convert(CuPtr{T}, buf.ptr)

function view(buf::Buffer, bytes::Int)
    bytes > buf.bytesize && throw(BoundsError(buf, bytes))
    return Mem.Buffer(buf.ptr+bytes, buf.bytesize-bytes, buf.ctx)
end

## refcounting

const refcounts = Dict{Buffer, Int}()

function refcount(buf::Buffer)
    get(refcounts, Base.unsafe_convert(CuPtr{Cvoid}, buf), 0)
end

"""
    retain(buf)

Increase the refcount of a buffer.
"""
function retain(buf::Buffer)
    refcount = get!(refcounts, buf, 0)
    refcounts[buf] = refcount + 1
    return
end

"""
    release(buf)

Decrease the refcount of a buffer. Returns `true` if the refcount has dropped to 0, and
some action needs to be taken.
"""
function release(buf::Buffer)
    haskey(refcounts, buf) || error("Release of unmanaged $buf")
    refcount = refcounts[buf]
    @assert refcount > 0 "Release of dead $buf"
    refcounts[buf] = refcount - 1
    return refcount==1
end

## memory info

"""
    info()

Returns a tuple of two integers, indicating respectively the free and total amount of memory
(in bytes) available for allocation by the CUDA context.
"""
function info()
    free_ref = Ref{Csize_t}()
    total_ref = Ref{Csize_t}()
    @apicall(:cuMemGetInfo, (Ptr{Csize_t},Ptr{Csize_t}), free_ref, total_ref)
    return convert(Int, free_ref[]), convert(Int, total_ref[])
end

"""
    free()

Returns the free amount of memory (in bytes), available for allocation by the CUDA context.
"""
free() = info()[1]

"""
    total()

Returns the total amount of memory (in bytes), available for allocation by the CUDA context.
"""
total() = info()[2]

"""
    used()

Returns the used amount of memory (in bytes), allocated by the CUDA context.
"""
used() = total()-free()

## generic interface (for documentation purposes)

"""
Allocate linear memory on the device and return a buffer to the allocated memory. The
allocated memory is suitably aligned for any kind of variable. The memory will not be freed
automatically, use [`free(::Buffer)`](@ref) for that.
"""
function alloc end

"""
Free device memory.
"""
function free end

"""
Initialize device memory with a repeating value.
"""
function set! end

"""
Upload memory from host to device.
Executed asynchronously on `stream` if `async` is true.
"""
function upload end
@doc (@doc upload) upload!

"""
Download memory from device to host.
Executed asynchronously on `stream` if `async` is true.
"""
function download end
@doc (@doc download) download!

"""
Transfer memory from device to device.
Executed asynchronously on `stream` if `async` is true.
"""
function transfer end
@doc (@doc transfer) transfer!

## pointer-based

@enum(CUmem_attach, ATTACH_GLOBAL = 0x01,
                    ATTACH_HOST   = 0x02)
                    #ATTACH_SINGLE = 0x04) # Defined but not valid

"""
    alloc(bytes::Integer)

Allocate `bytesize` bytes of memory.

Note that, contrary to the CUDA API, zero-size allocations are permitted. Such allocations
will point to the null pointer, and are not attached to a valid context.
"""
function alloc(bytesize::Integer, managed=false; flags::CUmem_attach=ATTACH_GLOBAL)
    bytesize == 0 && return Buffer(CU_NULL, 0, CuContext(C_NULL))

    ptr_ref = Ref{CuPtr{Cvoid}}()
    if !managed
        @apicall(:cuMemAlloc,
                 (Ptr{CuPtr{Cvoid}}, Csize_t),
                 ptr_ref, bytesize)
    else
        @apicall(:cuMemAllocManaged,
                 (Ptr{CuPtr{Cvoid}}, Csize_t, Cuint),
                 ptr_ref, bytesize, flags)
    end
    return Buffer(ptr_ref[], bytesize, CuCurrentContext())
end

function prefetch(buf::Buffer, bytes=buf.bytesize; stream::CuStream=CuDefaultStream())
    bytes > buf.bytesize && throw(BoundsError(buf, bytes))
    dev = device(buf.ctx)
    @apicall(:cuMemPrefetchAsync, (CuPtr{Cvoid}, Csize_t, CuDevice_t, CuStream_t),
             buf, bytes, dev, stream)
end

@enum(CUmem_advise, ADVISE_SET_READ_MOSTLY          = 0x01,
                    ADVISE_UNSET_READ_MOSTLY        = 0x02,
                    ADVISE_SET_PREFERRED_LOCATION   = 0x03,
                    ADVISE_UNSET_PREFERRED_LOCATION = 0x04,
                    ADVISE_SET_ACCESSED_BY          = 0x05,
                    ADVISE_UNSET_ACCESSED_BY        = 0x06)

function advise(buf::Buffer, advice::CUmem_advise, bytes=buf.bytesize, device=device(buf.ctx))
    bytes > buf.bytesize && throw(BoundsError(buf, bytes))
    @apicall(:cuMemAdvise, (CuPtr{Cvoid}, Csize_t, Cuint, CuDevice_t),
             buf, bytes, advice, device)
end

function free(buf::Buffer)
    if buf.ptr != C_NULL
        @apicall(:cuMemFree, (CuPtr{Cvoid},), buf)
    end
    return
end

for T in [UInt8, UInt16, UInt32]
    bits = 8*sizeof(T)
    fn_sync = Symbol("cuMemsetD$(bits)")
    fn_async = Symbol("cuMemsetD$(bits)Async")
    @eval begin
        @doc $"""
            set!(buf::Buffer, value::$T, len::Integer, [stream=CuDefaultStream()]; async=false)

        Initialize device memory by copying the $bits-bit value `val` for `len` times.
        Executed asynchronously if `async` is true.
        """
        function set!(buf::Buffer, value::$T, len::Integer,
                      stream::CuStream=CuDefaultStream(); async::Bool=false)
            if async
                @apicall($(QuoteNode(fn_async)),
                         (CuPtr{Cvoid}, $T, Csize_t, CuStream_t),
                         buf.ptr, value, len, stream)
            else
                @assert stream==CuDefaultStream()
                @apicall($(QuoteNode(fn_sync)),
                         (CuPtr{Cvoid}, $T, Csize_t),
                         buf.ptr, value, len)
            end
        end
    end
end

"""
    upload!(dst::Buffer, src, nbytes::Integer, [stream=CuDefaultStream()]; async=false)

Upload `nbytes` memory from `src` at the host to `dst` on the device.
"""
function upload!(dst::Buffer, src::Ref, nbytes::Integer,
                 stream::CuStream=CuDefaultStream(); async::Bool=false)
    if async
        @apicall(:cuMemcpyHtoDAsync,
                 (CuPtr{Cvoid}, Ptr{Cvoid}, Csize_t, CuStream_t),
                 dst, src, nbytes, stream)
    else
        @assert stream==CuDefaultStream()
        @apicall(:cuMemcpyHtoD,
                 (CuPtr{Cvoid}, Ptr{Cvoid}, Csize_t),
                 dst, src, nbytes)
    end
end

"""
    download!(dst::Ref, src::Buffer, nbytes::Integer, [stream=CuDefaultStream()]; async=false)

Download `nbytes` memory from `src` on the device to `src` on the host.
"""
function download!(dst::Ref, src::Buffer, nbytes::Integer,
                   stream::CuStream=CuDefaultStream(); async::Bool=false)
    if async
        @apicall(:cuMemcpyDtoHAsync,
                 (Ptr{Cvoid}, CuPtr{Cvoid}, Csize_t, CuStream_t),
                 dst, src, nbytes, stream)
    else
        @assert stream==CuDefaultStream()
        @apicall(:cuMemcpyDtoH,
                 (Ptr{Cvoid}, CuPtr{Cvoid}, Csize_t),
                 dst, src, nbytes)
    end
end

"""
    transfer!(dst::Buffer, src::Buffer, nbytes::Integer, [stream=CuDefaultStream()]; async=false)

Transfer `nbytes` of device memory from `src` to `dst`.
"""
function transfer!(dst::Buffer, src::Buffer, nbytes::Integer,
                   stream::CuStream=CuDefaultStream(); async::Bool=false)
    if async
        @apicall(:cuMemcpyDtoDAsync,
                 (CuPtr{Cvoid}, CuPtr{Cvoid}, Csize_t, CuStream_t),
                 dst, src, nbytes, stream)
    else
        @assert stream==CuDefaultStream()
        @apicall(:cuMemcpyDtoD,
                 (CuPtr{Cvoid}, CuPtr{Cvoid}, Csize_t),
                 dst, src, nbytes)
    end
end

## array based

"""
    alloc(src::AbstractArray)

Allocate space to store the contents of `src`.
"""
function alloc(src::AbstractArray)
    return alloc(sizeof(src))
end

"""
    upload!(dst::Buffer, src::AbstractArray, [stream=CuDefaultStream()]; async=false)

Upload the contents of an array `src` to `dst`.
"""
function upload!(dst::Buffer, src::AbstractArray,
                 stream=CuDefaultStream(); async::Bool=false)
    upload!(dst, Ref(src, 1), sizeof(src), stream; async=async)
end

"""
    upload(src::AbstractArray)::Buffer

Allocates space for and uploads the contents of an array `src`, returning a Buffer.
Cannot be executed asynchronously due to the synchronous allocation.
"""
function upload(src::AbstractArray) # TODO: stream, async
    dst = alloc(src)
    upload!(dst, src)
    return dst
end

"""
    download!(dst::AbstractArray, src::Buffer, [stream=CuDefaultStream()]; async=false)

Downloads memory from `src` to the array at `dst`. The amount of memory downloaded is
determined by calling `sizeof` on the array, so it needs to be properly preallocated.
"""
function download!(dst::AbstractArray, src::Buffer,
                   stream::CuStream=CuDefaultStream(); async::Bool=false)
    ref = Ref(dst, 1)
    download!(ref, src, sizeof(dst), stream; async=async)
    return
end

## type based

function check_type(::Type{Buffer}, T)
    if isa(T, UnionAll) || T.abstract || !isconcretetype(T)
        throw(ArgumentError("cannot represent abstract or non-leaf object"))
    end
    Base.datatype_pointerfree(T) || throw(ArgumentError("cannot handle non-ptrfree objects"))
    sizeof(T) == 0 && throw(ArgumentError("cannot represent singleton objects"))
end

"""
    alloc(T::Type, [count::Integer=1])

Allocate space for `count` objects of type `T`.
"""
function alloc(::Type{T}, count::Integer=1) where {T}
    check_type(Buffer, T)

    return alloc(sizeof(T)*count)
end

"""
    download(::Type{T}, src::Buffer, [count::Integer=1], [stream=CuDefaultStream()]; async=false)::Vector{T}

Download `count` objects of type `T` from the device at `src`, returning a vector.
"""
function download(::Type{T}, src::Buffer, count::Integer=1,
                  stream::CuStream=CuDefaultStream(); async::Bool=false) where {T}
    dst = Vector{T}(undef, count)
    download!(dst, src, stream; async=async)
    return dst
end

end

ext: jl


text:
module MEstimation # general estimating equations with or without bias-reducting adjustments

using NLsolve
using Optim
using FiniteDiff
using ForwardDiff
using LinearAlgebra
using Distributions
using ReverseDiff
# using InvertedIndices

import Base: show, print
import StatsBase: fit, aic, vcov, coef, coeftable, stderror, CoefTable

export objective_function
export objective_function_template

export estimating_function
export get_estimating_function
export estimating_function_template

export aic
export tic
export vcov
export coef
export coeftable
export fit
export stderror

export slice
# export profile

const JCACHE = Dict{DataType,Any}()
const GCACHE = Dict{DataType,Any}()
const EOCACHE = Dict{DataType,Any}()
# look at objective_functions.jl
# take a test and try to rewrite the test using the new design (look for all the new things I will need to do)
# remove most of the tests, keep only logreg using objectives
include("estimating_functions.jl")
include("objective_functions.jl")
include("fit.jl")
include("result_methods.jl")
include("slice.jl")
# include("profile.jl")

end # module


text:
module MetroMapRuby
    using JuMP
    import MathOptInterface
    using DataStructures
    import MathProgBase

    include("types.jl")
    include("transitmap.jl")
    include("mipmodel.jl")
    include("optimze.jl")

    export Station
    export Edge
    export GeoCoordinate
    export GeoStation
    export EuclideanCoordinate
    export EuclideanStation
    export EuclideanEdge
    export Line
    export GeoTransitMap
    export InputGraph
    export ProcessedEdge
    export optimize
    export TransitMapLayout
    export angle_deg
    export reduce_transitmap
    export classify_direction_sector
    export restore_transitmap
end # module

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
using FunctionalCollections
using Test

@testset "Persistent Queues" begin

    @testset "equality" begin
        @test PersistentQueue(1:100) == PersistentQueue(1:100)
        @test PersistentQueue(1:100) != (PersistentQueue(1:99))
    end

    @testset "isempty" begin
        @test isempty(PersistentQueue{Int}())
        @test !isempty(PersistentQueue([1]))
    end

    @testset "peek" begin
        @test peek(PersistentQueue(1:100)) == 100
        @test try peek(PersistentQueue{Int}()); false catch e true end
    end

    @testset "pop" begin
        @test pop(PersistentQueue(1:100)) == PersistentQueue(1:99)
        @test try pop(PersistentQueue{Int}()); false catch e true end
    end

    @testset "enq" begin
        q = PersistentQueue{Int}()
        @test peek(enq(q, 1)) == 1
        @test peek(pop(pop(enq(enq(enq(q, 1), 2), 3)))) == 3
    end

    @testset "iteration" begin
        arr2 = Int[]
        for i in PersistentQueue(1:1000)
            push!(arr2, i)
        end
        @test arr2 == collect(1000:-1:1)
    end

    @testset "hash" begin
        @test hash(PersistentQueue(1:1000)) == hash(PersistentQueue(1:1000))
    end

    @testset "length" begin
        @test length(PersistentQueue([1, 2, 3])) == 3
        @test length(PersistentQueue{Int}()) == 0

        @test length(pop(PersistentQueue([1, 2, 3]))) == 2
    end

end


text:
function AngleToX360(v)
    θ = angle(Complex(v...))
    if θ < 0
        θ = 2*π+θ
    end
    return θ
end

function Angle(u,v)
    α = (u⋅v)/(norm(u)*norm(v))
    if abs(α) > 1
        α=sign(α)*1
    end
    return acos(α)
end

function unwrap(v, inplace=false)
  # currently assuming an array
  unwrapped = inplace ? v : copy(v)
  for i in 2:length(v)
    while unwrapped[i] - unwrapped[i-1] >= pi/2
      unwrapped[i] -= pi
    end
    while unwrapped[i] - unwrapped[i-1] <= -pi/2
      unwrapped[i] += pi
    end
  end
  return unwrapped
end

function Ellipse(n;a=1,b=1,max=2π,θ=0,center=[0.0,0.0])
    x = [a*cos(t)*cos(θ) - b*sin(t)*sin(θ) for t in 0:(2π/n):2π] .+ center[1]
    y = [a*cos(t)*sin(θ) + b*sin(t)*cos(θ) for t in 0:(2π/n):2π] .+ center[2]
    return cat(x,y,dims=2)
end

beetle_length = 13.0
beetle_width = 7.0
dt = 1.0/30.0

function Density(X)
    tess,inds = AlphaShapes.GetDelaunayTriangulation(X,true);
	d = zeros(size(X,1))
	for i in 1:size(X,1)
    	d[i] = AlphaShapes.WeightedDTFELocalDensity(i,tess,inds)
	end
	return d
end

flatten(x) = reshape(x,prod(size(x)))

function logspace(start,stop,n=10)
    s = log10(start)
    S = log10(stop)
    x = []
    for i in linspace(s,S,n)
        push!(x,10^i)
    end
    return x
end

function linspace(l,u,n=10)
    return collect(range(l,stop=u,length=n))
end

function DVBinPlot(;D,V,n_bins=10,beetlescale=true,beetle_length=14.0,xscale=:log10,FPS=30.0,errors=true)
    IFT = 1/FPS
    bin_values = logspace(minimum(D),maximum(D),n_bins)

    Bin_V_SE = [[] for i in 1:n_bins]
    Bin_V = zeros(n_bins)
    number = zeros(n_bins)

    for b in 1:(n_bins-1)
        ind1 = findall(x->x.>=bin_values[b],D)
        v = V[ind1]
        ind = findall(x->x.<bin_values[b+1],D[ind1])
        if size(ind,1) > 0
            Bin_V[b] += sum(v[ind])
            append!(Bin_V_SE[b],v[ind])
            number[b] += size(ind,1)
        end
    end

    if beetlescale
        SE = zeros(n_bins)
        for i in 1:n_bins
            if size(Bin_V_SE[i],1) > 0
                SE[i] = std(Float64.(Bin_V_SE[i]) ./ (beetle_length * IFT)) / sqrt(number[i])
            end
        end
        Bin_V = Bin_V[2:end-1] ./ ((number[2:end-1] .* beetle_length) * IFT)
        bin_values = bin_values[2:end-1].*beetle_length^2.0
        if errors
            scatter(bin_values,Bin_V,xscale=xscale,yerr=SE,label="")
        else
            plot(bin_values,Bin_V,xscale=xscale,label="")
        end
    else
        SE = zeros(n_bins)
        for i in 1:n_bins
            if size(Bin_V_SE[i],1) > 0
                SE[i] = std(Float64.(Bin_V_SE[i])) / sqrt(number[i])
            end
        end
        if errors
            plot(bin_values,Bin_V,xscale=xscale,yerr=SE,label="")
        else
            plot(bin_values,Bin_V,xscale=xscale,label="")
        end
    end
    println(number)
    return bin_values,Bin_V,SE[2:end-1]
end

function MSD(X,t=size(X,1),lags=120)
    l = 0:lags
    D = zeros(size(l,1))
    C = zeros(size(l,1),size(X,1))
    @showprogress for i in 1:t
            for τ in l
                if i+τ > size(X,1)
                    continue
                end
                for j in 1:size(X,2)
                    p = norm(X[i+τ,j,1:2] - X[i,j,1:2])^2
                    if isnan(p) == false
                        D[τ+1] += p
                        C[τ+1,i] += 1
                    end
                end
            end
        end
    return D ./ sum(C,dims=2),collect(l)
end

function msdθ(X,t=size(X,1),lags=30)
    l = 0:lags
    D = zeros(size(l,1))
    A = X[:,:,3]
    C = zeros(size(l,1),t)
    @showprogress for i in 1:t
            for τ in l
                if i+τ > size(A,1)
                    continue
                end
                for j in 1:size(X,2)
                    p = (A[i+τ,j] - A[i,j])^2
                    if isnan(p) == false
                        D[τ+1] += p
                        C[τ+1,i] += 1
                    end
                end
            end
        end
    return D ./ sum(C,dims=2),collect(l)
end

function VelocityHeading(X,t,b)
    h = X[t+1,b,1:2].-X[t,b,1:2]
    return h./norm(h)
end

function Arrow(d,pos,K=1,c="red",linewidth=1)
    v = ones(2).*d*K + pos
    plot!([pos[1],v[1]],[pos[2],v[2]],color=c,label="",linewidth=linewidth,arrow=arrow())
end

function TimeBinDensity(BX,bins,T;length_scale=2.0)
    times = Int.(floor.(range(1,stop=size(T,1),length=bins)))
    collect(times)
    pdfs = zeros(bins-1,length(BX))
    for t in 1:size(times,1)-1
        start = times[t]
        stop = times[t+1]
        Y = T[start:stop,:,1:2]

        D = zeros(size(Y,1),size(Y,2))
        for i in 1:size(Y,1)
            D[i,:] = Density(Y[i,:,1:2]./length_scale)
        end
        D = replace(D,NaN => 0.0)
        x = flatten(D)
        x = x[x.>0]
        x = x[x.<Inf]

        h = fit(Histogram,x,BX)
        h = normalize(h, mode=:density)
        h.weights = h.weights./size(x,1)
        push!(h.weights,0.0)
        pdfs[t,:] = h.weights
    end
    return pdfs
end

function TimeBinDistribution(x,data;bins=10)
    times = Int.(floor.(range(1,stop=size(data,1),length=bins)))
    collect(times)
    pdfs = zeros(bins-1,length(x))
    for t in 1:size(times,1)-1
        start = times[t]
        stop = times[t+1]
        Y = []
        for i in start:stop
            for j in 1:size(data[i],1)
                push!(Y,data[i][j])
            end
        end
        y = replace!(Y,NaN => 0.0)
        y = y[y.>0]
        y = y[y.<Inf]

		h = fit(Histogram,y,x)
		h = normalize(h, mode=:density)
		h.weights = h.weights./size(x,1)
		push!(h.weights,0.0)
		pdfs[t,:] = h.weights
    end
    return pdfs
end

## Collision Filtering

using Distances

function Collisions(R=13.0,data=[X200])
    R = 13.0
    Data = []
    Cs = []
    for X in data
        C = Bool.(zeros(size(X,1),size(X,2),size(X,2)))
        for t in 1:size(X,1)
            d = pairwise(Euclidean(),X[t,:,1:2]').<=R
            for b in 1:size(X,2)
                for j in findall(x->x.==1,d[b,:])
                    C[t,b,j] = 1
                end
                C[t,b,b] = 0
            end
        end

        push!(Cs,C)

        for b in 1:size(C,2)
            flag = sum(C[1,b,:]) > 0
            d = []
            for t in 2:size(C,1)
                if sum(C[t,b,:]) > 0
                    push!(d,[t,X[t,b,:]...,findall(x->x.==1,C[t,b,:])...])
                    flag = 0
                end
                if sum(C[t,b,:]) == 0
                    if flag == 0
                        push!(Data,d)
                        d = []
                    end
                    flag = 1
                end
            end
            push!(Data,d)
        end
    end
    @info length(Data)
    return Cs
end

function CollisionFreeTrajectories(C,X)
    T = []
    d = zeros(size(X,1),size(X,2))
    for i in 1:size(X,1)
        d[i,:] = Density(X[i,:,1:2]./13.0)
    end
    @showprogress for i in 1:size(C,2)
        c = sum(C[:,i,:],dims=2)[:,1]
        flag = c[1]
        for j in 1:size(c,1)
            tau = []
            offset = 0
            while length(findall(x->x.==0,c)) > 0
                s = findfirst(x->x.==0,c)
                e = findfirst(x->x.==1,c)
                if e == nothing
                    e = length(c)
                end
                for k in s:1:e
                    push!(tau,([X[offset+k,i,:]...,d[offset+k,i]],offset+k))
                end
                offset += e
                c = c[(e+1):end]
            end
            if length(tau) > 0
                push!(T,tau)
            end
        end
    end

    # gets length count distribution
    L = []
    O = []
    for i in 1:size(T,1)
        t = T[i]
        if length(t) > 1
            o = [t[1][1]]
            k = t[1][2]
            l = 0
            for i in 2:size(t,1)
                if t[i][2] - k == 1
                    l += 1
                    k += 1
                    push!(o,t[i][1])
                else
                    k = t[i][2]
                    push!(L,l)
                    l = 0
                    o = [t[i][1]]
                    push!(O,o)
                end
            end
        end
    end
    return T,L,O
end

function MSDθ_CF(data,dt=1.0/30.0,lag_max=30)
    C = [[] for i in 1:(lag_max+1)]
    @showprogress for i in 1:length(data)
        d = data[i]
        if length(d) >= 2
            for t in 1:size(d,1)
                for (i,tau) in enumerate(0:lag_max)
                    if 1 < t+tau < size(d,1)
                        c = (d[t+tau][3] - d[t][3])^2.0
                        if isnan(c) == false
                            push!(C[i],c)
                        end
                    end
                end
            end
        end
    end
    return C,collect(0:lag_max).*dt
end

function MSD_CF(data,dt=1.0/30.0,lag_max=30)
    C = [[] for i in 1:(lag_max+1)]
    @showprogress for i in 1:length(data)
        d = data[i]
        if length(d) >= 2
            for t in 1:size(d,1)
                for (i,tau) in enumerate(0:lag_max)
                    if 1 < t+tau < size(d,1)
                        c = (norm(d[t+tau][1:2] .- d[t][1:2]))^2.0
                        if isnan(c) == false
                            push!(C[i],c)
                        end
                    end
                end
            end
        end
    end
    return C,collect(0:lag_max).*dt
end

function NVCorr(data,dt=1.0/30.0,lag_max=30,min_length=2,info=false)
    C = [[] for i in 1:(2*lag_max+1)]
    n = 0
    @showprogress for i in 1:length(data)
        d = data[i]
        if length(d) >= min_length
            n += 1
			# central difference
            for t in 1:size(d,1)
                for (i,tau) in enumerate(-lag_max:1:lag_max)
                    if 1 < t+tau < size(d,1)-1
						# central difference about t+tau
                        h = d[t+tau+1][1:2] .- d[t+tau-1][1:2]
                        h = h ./ norm(h)
                        o = [cos(d[t][3]),sin(d[t][3])]
                        c = h⋅o
                        if isnan(c) == false
                            push!(C[i],c)
                        end
                    end
                end
            end
        end
    end
    if info
        @info n
    end
    return C,collect(-lag_max:1:lag_max).*dt
end

max_stars_repo_name: harveydevereux/CUDA-Whirligigs

max_issues_repo_name: harveydevereux/CUDA-Whirligigs

max_forks_repo_name: harveydevereux/CUDA-Whirligigs


text:
using DashBootstrapComponents

popovers = html_div([
    # First example - using dbc_popoverbody
    dbc_button("popover-target", id = "popover-target", className = "me-1"),
    dbc_popover(
        dbc_popoverbody("My `target` is `popover-target`."),
        target = "popover-target",
        trigger = "click",
    ),
    # Second example - using body=true
    dbc_button("popover-alt-target", id = "popover-alt-target", className = "me-1"),
    dbc_popover(
        "My `target` is `popover-alt-target`.",
        body = true,
        target = "popover-alt-target",
        trigger = "click",
    ),
])

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
### This support BLAS style multiplication
#           α * A * B + β C
# but avoids the broadcast machinery

# Lazy representation of α*A*B + β*C
struct MulAdd{StyleA, StyleB, StyleC, T, AA, BB, CC}
    α::T
    A::AA
    B::BB
    β::T
    C::CC
end

@inline MulAdd{StyleA,StyleB,StyleC}(α::T, A::AA, B::BB, β::T, C::CC) where {StyleA,StyleB,StyleC,T,AA,BB,CC} =
    MulAdd{StyleA,StyleB,StyleC,T,AA,BB,CC}(α,A,B,β,C)

@inline function MulAdd{StyleA,StyleB,StyleC}(αT, A, B, βV, C) where {StyleA,StyleB,StyleC}
    α,β = promote(αT,βV)
    MulAdd{StyleA,StyleB,StyleC}(α, A, B, β, C)
end

@inline MulAdd(α, A::AA, B::BB, β, C::CC) where {AA,BB,CC} =
    MulAdd{typeof(MemoryLayout(AA)), typeof(MemoryLayout(BB)), typeof(MemoryLayout(CC))}(α, A, B, β, C)

MulAdd(A, B) = MulAdd(Mul(A, B))
function MulAdd(M::Mul)
    TV = eltype(M)
    MulAdd(scalarone(TV), M.A, M.B, scalarzero(TV), mulzeros(TV,M))
end

@inline eltype(::MulAdd{StyleA,StyleB,StyleC,T,AA,BB,CC}) where {StyleA,StyleB,StyleC,T,AA,BB,CC} =
     promote_type(_mul_eltype(T, eltype(AA), eltype(BB)), _mul_eltype(T, eltype(CC)))

size(M::MulAdd, p::Int) = size(M)[p]
axes(M::MulAdd, p::Int) = axes(M)[p]
length(M::MulAdd) = prod(size(M))
size(M::MulAdd) = map(length,axes(M))
axes(M::MulAdd) = axes(M.C)

similar(M::MulAdd, ::Type{T}, axes) where {T,N} = similar(Array{T}, axes)
similar(M::MulAdd{<:DualLayout,<:Any,<:Any,<:Any,<:Adjoint}, ::Type{T}, axes) where {T,N} = similar(Array{T}, axes[2])'
similar(M::MulAdd{<:DualLayout,<:Any,<:Any,<:Any,<:Transpose}, ::Type{T}, axes) where {T,N} = transpose(similar(Array{T}, axes[2]))
similar(M::MulAdd, ::Type{T}) where T = similar(M, T, axes(M))
similar(M::MulAdd) = similar(M, eltype(M))

function checkdimensions(M::MulAdd)
    @boundscheck check_mul_axes(M.α, M.A, M.B)
    @boundscheck check_mul_axes(M.β, M.C)
    @boundscheck axes(M.A,1) == axes(M.C,1) || throw(DimensionMismatch("First axis of A, $(axes(M.A,1)), and first axis of C, $(axes(M.C,1)) must match"))
    @boundscheck axes(M.B,2) == axes(M.C,2) || throw(DimensionMismatch("Second axis of B, $(axes(M.B,2)), and second axis of C, $(axes(M.C,2)) must match"))
end
@propagate_inbounds function instantiate(M::MulAdd)
    checkdimensions(M)
    M
end

const ArrayMulArrayAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractArray,<:AbstractArray,<:AbstractArray}
const MatMulVecAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractMatrix,<:AbstractVector,<:AbstractVector}
const MatMulMatAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractMatrix,<:AbstractMatrix,<:AbstractMatrix}
const VecMulMatAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractVector,<:AbstractMatrix,<:AbstractMatrix}

broadcastable(M::MulAdd) = M

const BlasMatMulVecAdd{StyleA,StyleB,StyleC,T<:BlasFloat} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractMatrix{T},<:AbstractVector{T},<:AbstractVector{T}}
const BlasMatMulMatAdd{StyleA,StyleB,StyleC,T<:BlasFloat} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractMatrix{T},<:AbstractMatrix{T},<:AbstractMatrix{T}}
const BlasVecMulMatAdd{StyleA,StyleB,StyleC,T<:BlasFloat} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractVector{T},<:AbstractMatrix{T},<:AbstractMatrix{T}}

muladd!(α, A, B, β, C) = materialize!(MulAdd(α, A, B, β, C))
materialize(M::MulAdd) = copy(instantiate(M))
copy(M::MulAdd) = copyto!(similar(M), M)

_fill_copyto!(dest, C) = copyto!(dest, C)
_fill_copyto!(dest, C::Zeros) = zero!(dest) # exploit special fill! overload

@inline copyto!(dest::AbstractArray{T}, M::MulAdd) where T =
    muladd!(M.α, unalias(dest,M.A), unalias(dest,M.B), M.β, _fill_copyto!(dest, M.C))

# Modified from LinearAlgebra._generic_matmatmul!
function tile_size(T, S, R)
    tile_size = 0
    if isbitstype(R) && isbitstype(T) && isbitstype(S)
        tile_size = floor(Int, sqrt(tilebufsize / max(sizeof(R), sizeof(S), sizeof(T))))
    end
    tile_size
end

function tiled_blasmul!(tile_size, α, A::AbstractMatrix{T}, B::AbstractMatrix{S}, β, C::AbstractMatrix{R}) where {S,T,R}
    mA, nA = size(A)
    mB, nB = size(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    size(C) == (mA, nB) || throw(DimensionMismatch("Dimensions must match"))

    (iszero(mA) || iszero(nB)) && return C
    iszero(nA) && return lmul!(β, C)

    @inbounds begin
        sz = (tile_size, tile_size)
        # FIXME: This code is completely invalid!!!
        Atile = Array{T}(undef, sz)
        Btile = Array{S}(undef, sz)

        z1 = zero(A[1, 1]*B[1, 1] + A[1, 1]*B[1, 1])
        z = convert(promote_type(typeof(z1), R), z1)

        if mA < tile_size && nA < tile_size && nB < tile_size
            copy_transpose!(Atile, 1:nA, 1:mA, 'N', A, 1:mA, 1:nA)
            copyto!(Btile, 1:mB, 1:nB, 'N', B, 1:mB, 1:nB)
            for j = 1:nB
                boff = (j-1)*tile_size
                for i = 1:mA
                    aoff = (i-1)*tile_size
                    s = z
                    for k = 1:nA
                        s += Atile[aoff+k] * Btile[boff+k]
                    end
                    C[i,j] = α*s + β*C[i,j]
                end
            end
        else
            Ctile = Array{R}(undef, sz)
            for jb = 1:tile_size:nB
                jlim = min(jb+tile_size-1,nB)
                jlen = jlim-jb+1
                for ib = 1:tile_size:mA
                    ilim = min(ib+tile_size-1,mA)
                    ilen = ilim-ib+1
                    copyto!(Ctile, 1:ilen, 1:jlen, C, ib:ilim, jb:jlim)
                    lmul!(β,Ctile)
                    for kb = 1:tile_size:nA
                        klim = min(kb+tile_size-1,mB)
                        klen = klim-kb+1
                        copy_transpose!(Atile, 1:klen, 1:ilen, 'N', A, ib:ilim, kb:klim)
                        copyto!(Btile, 1:klen, 1:jlen, 'N', B, kb:klim, jb:jlim)
                        for j=1:jlen
                            bcoff = (j-1)*tile_size
                            for i = 1:ilen
                                aoff = (i-1)*tile_size
                                s = z
                                for k = 1:klen
                                    s += Atile[aoff+k] * Btile[bcoff+k]
                                end
                                Ctile[bcoff+i] += α*s
                            end
                        end
                    end
                    copyto!(C, ib:ilim, jb:jlim, Ctile, 1:ilen, 1:jlen)
                end
            end
        end
    end

    C
end

function default_blasmul!(α, A::AbstractMatrix, B::AbstractMatrix, β, C::AbstractMatrix)
    mA, nA = size(A)
    mB, nB = size(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    size(C) == (mA, nB) || throw(DimensionMismatch("Dimensions must match"))

    lmul!(β, C)

    (iszero(mA) || iszero(nB)) && return C
    iszero(nA) && return C

    @inbounds for k in colsupport(A), j in rowsupport(B,rowsupport(A,k))
        z2 = zero(A[k, 1]*B[1, j] + A[k, 1]*B[1, j])
        Ctmp = convert(promote_type(eltype(C), typeof(z2)), z2)
        @simd for ν = rowsupport(A,k) ∩ colsupport(B,j)
            Ctmp = muladd(A[k, ν],B[ν, j],Ctmp)
        end
        C[k,j] = muladd(α,Ctmp, C[k,j])
    end
    C
end

function _default_blasmul!(::IndexLinear, α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector)
    mA, nA = size(A)
    mB = length(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    length(C) == mA || throw(DimensionMismatch("Dimensions must match"))

    lmul!(β, C)
    (nA == 0 || mB == 0)  && return C

    z = zero(A[1]*B[1] + A[1]*B[1])
    Astride = size(A, 1) # use size, not stride, since its not pointer arithmetic

    @inbounds for k in colsupport(B,1)
        aoffs = (k-1)*Astride
        b = B[k]
        for i = 1:mA
            C[i] += α * A[aoffs + i] * b
        end
    end

    C
end

function _default_blasmul!(::IndexCartesian, α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector)
    mA, nA = size(A)
    mB = length(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    length(C) == mA || throw(DimensionMismatch("Dimensions must match"))

    lmul!(β, C)
    (nA == 0 || mB == 0)  && return C

    z = zero(A[1,1]*B[1] + A[1,1]*B[1])

    @inbounds for k in colsupport(B,1)
        b = B[k]
        for i = colsupport(A,k)
            C[i] += α * A[i,k] * b
        end
    end

    C
end

default_blasmul!(α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector) =
    _default_blasmul!(Base.IndexStyle(typeof(A)), α, A, B, β, C)

function materialize!(M::MatMulMatAdd)
    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
    default_blasmul!(α, unalias(C,A), unalias(C,B), iszero(β) ? false : β, C)
end

function materialize!(M::MatMulMatAdd{<:AbstractStridedLayout,<:AbstractStridedLayout,<:AbstractStridedLayout})
    α, Ain, Bin, β, C = M.α, M.A, M.B, M.β, M.C
    A = unalias(C, Ain)
    B = unalias(C, Bin)
    ts = tile_size(eltype(A), eltype(B), eltype(C))
    if iszero(β) # false is a "strong" zero to wipe out NaNs
        if ts == 0 || !(axes(A) isa NTuple{2,OneTo{Int}}) || !(axes(B) isa NTuple{2,OneTo{Int}}) || !(axes(C) isa NTuple{2,OneTo{Int}})
            default_blasmul!(α, A, B, false, C)
        else
            tiled_blasmul!(ts, α, A, B, false, C)
        end
    else
        if ts == 0 || !(axes(A) isa NTuple{2,OneTo{Int}}) || !(axes(B) isa NTuple{2,OneTo{Int}}) || !(axes(C) isa NTuple{2,OneTo{Int}})
            default_blasmul!(α, A, B, β, C)
        else
            tiled_blasmul!(ts, α, A, B, β, C)
        end
    end
end

function materialize!(M::MatMulVecAdd)
    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
    default_blasmul!(α, unalias(C,A), unalias(C,B), iszero(β) ? false : β, C)
end

@inline _gemv!(tA, α, A, x, β, y) = BLAS.gemv!(tA, α, unalias(y,A), unalias(y,x), β, y)
@inline _gemm!(tA, tB, α, A, B, β, C) = BLAS.gemm!(tA, tB, α, unalias(C,A), unalias(C,B), β, C)

# work around pointer issues
@inline materialize!(M::BlasMatMulVecAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractStridedLayout}) =
    _gemv!('N', M.α, M.A, M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulVecAdd{<:AbstractRowMajor,<:AbstractStridedLayout,<:AbstractStridedLayout}) =
    _gemv!('T', M.α, transpose(M.A), M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulVecAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _gemv!('C', M.α, adjoint(M.A), M.B, M.β, M.C)

@inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
    _gemm!('N', 'N', M.α, M.A, M.B, M.β, M.C)
@inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
    _gemm!('N', 'T', M.α, M.A, transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('N', 'C', M.α, M.A, adjoint(M.B), M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
    _gemm!('N', 'N', M.α, M.A, M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
    _gemm!('N', 'T', M.α, M.A, transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('N', 'C', M.α, M.A, adjoint(M.B), M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
    _gemm!('T', 'N', M.α, transpose(M.A), M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'N', M.α, adjoint(M.A), M.B, M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
    _gemm!('T', 'T', M.α, transpose(M.A), transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('T', 'C', M.α, transpose(M.A), adjoint(M.B), M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractRowMajor,<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'T', M.α, adjoint(M.A), transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'C', M.α, adjoint(M.A), adjoint(M.B), M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractRowMajor}) =
    _gemm!('T', 'T', M.α, M.B, M.A, M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('C', 'C', M.α, M.B, M.A, M.β, adjoint(M.C))

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractRowMajor}) =
    _gemm!('N', 'T', M.α, transpose(M.B), M.A, M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'T', M.α, transpose(M.B), M.A, M.β, M.C')
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'C', M.α, adjoint(M.B), M.A, M.β, adjoint(M.C))

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractColumnMajor,<:AbstractRowMajor}) =
    _gemm!('T', 'N', M.α, M.B, transpose(M.A), M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('C', 'N', M.α, M.B, adjoint(M.A), M.β, adjoint(M.C))

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractRowMajor,<:AbstractRowMajor}) =
    _gemm!('N', 'N', M.α, transpose(M.B), transpose(M.A), M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'N', M.α, adjoint(M.B), adjoint(M.A), M.β, adjoint(M.C))

###
# Symmetric
###

# make copy to make sure always works
@inline _symv!(tA, α, A, x, β, y) = BLAS.symv!(tA, α, unalias(y,A), unalias(y,x), β, y)
@inline _hemv!(tA, α, A, x, β, y) = BLAS.hemv!(tA, α, unalias(y,A), unalias(y,x), β, y)

materialize!(M::BlasMatMulVecAdd{<:SymmetricLayout{<:AbstractColumnMajor},<:AbstractStridedLayout,<:AbstractStridedLayout}) =
    _symv!(symmetricuplo(M.A), M.α, symmetricdata(M.A), M.B, M.β, M.C)

materialize!(M::BlasMatMulVecAdd{<:SymmetricLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout}) =
    _symv!(symmetricuplo(M.A) == 'L' ? 'U' : 'L', M.α, transpose(symmetricdata(M.A)), M.B, M.β, M.C)

materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:AbstractColumnMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _hemv!(symmetricuplo(M.A), M.α, hermitiandata(M.A), M.B, M.β, M.C)

materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _hemv!(symmetricuplo(M.A) == 'L' ? 'U' : 'L', M.α, hermitiandata(M.A)', M.B, M.β, M.C)

####
# Diagonal
####

# Diagonal multiplication never changes structure
similar(M::MulAdd{<:DiagonalLayout,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.B, T, axes)
similar(M::MulAdd{<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.B, T, axes)
similar(M::MulAdd{<:Any,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.A, T, axes)
# equivalent to rescaling
function materialize!(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}})
    checkdimensions(M)
    M.C .= (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C
    M.C
end

function materialize!(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}})
    checkdimensions(M)
    M.C .= M.α .* M.A .* getindex_value(M.B.diag) .+ M.β .* M.C
    M.C
end

BroadcastStyle(::Type{<:MulAdd}) = ApplyBroadcastStyle()

scalarone(::Type{T}) where T = one(T)
scalarone(::Type{<:AbstractArray{T}}) where T = scalarone(T)
scalarzero(::Type{T}) where T = zero(T)
scalarzero(::Type{<:AbstractArray{T}}) where T = scalarzero(T)

fillzeros(::Type{T}, ax) where T<:Number = Zeros{T}(ax)
mulzeros(::Type{T}, M) where T<:Number = fillzeros(T, axes(M))
mulzeros(::Type{T}, M::Mul{<:DualLayout,<:Any,<:Adjoint}) where T<:Number = fillzeros(T, axes(M,2))'
mulzeros(::Type{T}, M::Mul{<:DualLayout,<:Any,<:Transpose}) where T<:Number = transpose(fillzeros(T, axes(M,2)))

# initiate array-valued MulAdd
function _mulzeros!(dest::AbstractVector{T}, A, B) where T
    for k in axes(dest,1)
        dest[k] = similar(Mul(A[k,1],B[1]), eltype(T))
    end
    dest
end

function _mulzeros!(dest::AbstractMatrix{T}, A, B) where T
    for j in axes(dest,2), k in axes(dest,1)
        dest[k,j] = similar(Mul(A[k,1],B[1,j]), eltype(T))
    end
    dest
end

mulzeros(::Type{T}, M) where T<:AbstractArray = _mulzeros!(similar(Array{T}, axes(M)), M.A, M.B)

###
# Fill
###

copy(M::MulAdd{<:AbstractFillLayout,<:AbstractFillLayout,<:AbstractFillLayout}) = M.α*M.A*M.B + M.β*M.C

###
# DualLayout
###

transtype(::Adjoint) = adjoint
transtype(::Transpose) = transpose

function similar(M::MulAdd{<:DualLayout,<:Any,ZerosLayout}, ::Type{T}, (x,y)) where T
    @assert length(x) == 1
    trans = transtype(M.A)
    trans(similar(trans(M.A), T, y))
end

function similar(M::MulAdd{<:Any,<:DualLayout,ZerosLayout}, ::Type{T}, (x,y)) where T
    @assert length(x) == 1
    trans = transtype(M.B)
    trans(similar(trans(M.B), T, y))
end

const ZerosLayouts = Union{ZerosLayout,DualLayout{ZerosLayout}}
copy(M::MulAdd{<:ZerosLayouts, <:ZerosLayouts, <:ZerosLayouts}) = M.C
copy(M::MulAdd{<:ZerosLayouts, <:Any, <:ZerosLayouts}) = M.C
copy(M::MulAdd{<:Any, <:ZerosLayouts, <:ZerosLayouts}) = M.C


text:
"""
    @test_exists [ntests = 100] [maxtests = 1000] [logto = ""] argument_data, prop

Checks `prop` on randomly generated samples values specified in `argument_data`.
Returns `Pass` if `prop` holds for at least one value in the sample,
`Fail` if `prop` fails to hold for all sampled values,
and `Error` if a part of the code could not be evaluated
or if insufficient non-vacuous tests are generated.

## Description

`@test_exists` is to `@test_formany` as the existential quantifier
is to the universal quantifier.
That is,

    ( @test_exists x ∈ X, P(x)  ) == !( @test_formany x ∈ X, !P(x) )

Consequently, the same concerns with vacuity avoidance apply.
Additionally, analogous concerns with false positivity
apply: where `@test_formany` may falsely validate a property
by not exhausting its implicit test universe,
`@test_exists` may falsely *in*validate a property
by not exhausting its implicit search universe.

## Examples

```julia
julia> @test_exists ntests = 1000 -10 < x < 10, x^2 > 99
String["x = -9.966560160994264"]
Test Passed
  Expression: (:(ntests = 1000),:((-10 < x < 10,x ^ 2 > 99)))

## Here, the values of x, y, and z are generated independently,
## so a conditional statement with an antecedent that checks for
## order is necessary:
julia> @test_exists ntests = 1000 0 < x < 2*pi, 0 < y < 2*pi, 0 < z < 2*pi, x < y < z --> sin(x) < sin(y) < sin(z)
String["x = 0.3082312611033726","y = 0.49557010329840206","z = 2.5966119472766076"]
Test Passed
  Expression: (:((0 < x < 2pi,0 < y < 2pi,0 < z < 2pi,$(Expr(:-->, :(x < y < z), :(sin(x) < sin(y) < sin(z)))))),:(ntests = 1000))

## By contrast, here, x, y, and z are generated with sequential
## reference to the previously generated value(s),
## so they satisfy x<y<z automatically.
## No conditional is needed.
julia> @test_exists 0<x<2*pi, x<y<2*pi, y<z<2*pi, sin(x)<sin(y)<sin(z)
String["x = 5.787034350107664","y = 6.23560263220101","z = 6.26802436856299"]
Test Passed
  Expression: (:((0 < x < 2pi,x < y < 2pi,y < z < 2pi,sin(x) < sin(y) < sin(z))),)

## This returns Error, because the antecedent is too restrictive
## to generate 100 non-vacuous tests:
@test_exists ntests = 100 0 < x < 1000, x > 990 --> x > 999.9
Error During Test
  Test threw an exception of type ErrorException
  Expression: (:((0 < x < 1000,$(Expr(:-->, :(x > 990), :(x > 999.9))))),:(ntests = 100))
  Found only 14/100 values satisfying given condition.
```
"""
macro test_exists(exprs...)
    outex = Expr(:macrocall, Symbol("@test_cases"))
    for ex in exprs
        push!(outex.args, esc(ex))
    end
    ex1 = Expr(:(=), :mode, :test_exists)
    push!(outex.args, esc(ex1))
    return outex
end


text:
## File containing different utility functions ##

## Jittering object adapting to the type of the GP ##
struct Jittering end;

const jitter = Jittering()

@inline Base.Float64(::Jittering) = Float64(1e-5)
@inline Base.Float32(::Jittering) = Float32(1e-4)
@inline Base.Float16(::Jittering) = Float16(1e-3)
@inline Base.convert(::Type{Float64},::Jittering) = Float64(1e-5)
@inline Base.convert(::Type{Float32},::Jittering) = Float32(1e-4)
@inline Base.convert(::Type{Float16},::Jittering) = Float16(1e-3)

## delta function `(i,j)`, equal `1` if `i == j`, `0` else ##
@inline function δ(i::Integer,j::Integer)
    i == j ? 1.0 : 0.0
end

## Hadamard product between two arrays of same size ##
@inline function hadamard(A::AbstractArray{<:Real},B::AbstractArray{<:Real})
    A.*B
end

## Add on place the transpose of a matrix ##
@inline function add_transpose!(A::AbstractMatrix{<:Real})
    A .+= A'
end

## Return the trace of A*B' ##
@inline function opt_trace(A::AbstractMatrix{<:Real},B::AbstractMatrix{<:Real})
    dot(A,B)
end

## Return the diagonal of A*B' ##
@inline function opt_diag(A::AbstractArray{T,N},B::AbstractArray{T,N}) where {T<:Real,N}
    vec(sum(A.*B,dims=2))
end

## Return the multiplication of Diagonal(v)*B ##
function opt_diag_mul_mat(v::AbstractVector{T},B::AbstractMatrix{T}) where {T<:Real}
    v.*B
end

## Return the multiplication of κᵀ*Diagonal(θ)*κ
@inline function κdiagθκ(κ::AbstractMatrix{T},θ::AbstractVector{T}) where {T<:Real}
    transpose(θ.*κ)*κ
end

## Return the multiplication of ρ*κᵀ*Diagonal(θ)*κ
@inline function ρκdiagθκ(ρ::T,κ::AbstractMatrix{T},θ::AbstractVector{T}) where {T<:Real}
    transpose((ρ*θ).*κ)*κ
end

## Return the addition of a diagonal to a symmetric matrix ##
function opt_add_diag_mat(v::AbstractVector{T},B::AbstractMatrix{T}) where {T<:Real}
    A = copy(B)
    @inbounds for i in 1:size(A,1)
        A[i,i] += v[i]
    end
    A
end

## Temp fix until the deepcopy of the main package is fixed
function Base.copy(opt::Optimizer)
    f = length(fieldnames(typeof(opt)))
    copied_params = [deepcopy(getfield(opt, k)) for k = 1:f]
    return typeof(opt)(copied_params...)
end

Base.deepcopy(opt::Optimizer) = copy(opt)

## Compute exp(μ)/cosh(c) safely if there is an overflow ##
function safe_expcosh(μ::Real,c::Real)
    return isfinite(exp(μ)/cosh(c)) ? exp(μ)/cosh(c) : 2*logistic(2.0*max(μ,c))
end

## Return a safe version of log(cosh(c)) ##
function logcosh(c::Real)
    return log(exp(-2.0*c)+1.0)+c-logtwo
end

function logisticsoftmax(f::AbstractVector{<:Real})
    return normalize!(logistic.(f),1)
end

function logisticsoftmax(f::AbstractVector{<:Real},i::Integer)
    return logisticsoftmax(f)[i]
end

export make_grid
function make_grid(range1,range2)
    return hcat([i for i in range1, j in range2][:],[j for i in range1, j in range2][:])
end

function vech(A::AbstractMatrix{T}) where T
    m = size(A,1)
    v = Vector{T}(undef,m*(m+1)÷2)
    k = 0
    for j = 1:m, i = j:m
        @inbounds v[k += 1] = A[i,j]
    end
    return v
end

function duplicate_matrix(n::Int)
    m   = (n * (n + 1)) ÷ 2
    nsq = n^2
    D = zeros(Float64,nsq, m)
    row = 1
    a   = 1
    for i in 1:n
        b = i
        for j = 0:i-2
            D[row + j, b] = 1;
            b = b + n - j - 1;
        end
       row = row + i - 1

       for j = 0:n-i
           D[row + j, a + j] = 1
       end
       row = row + n - i + 1
       a   = a + n - i + 1
    end
    return D
end

function eliminate_matrix(m::Int)
    T = tril(ones(m,m))
    f = findall(T[:].==true)
    k = (m*(m+1))÷2
    m2 = m*m
    L = falses(k,m2)
    x = vec(f .+ m2*(0:(k-1)))
    L[x] .= true
    copy(transpose(L))
end

function commute_transpose(n::Int,m::Int=n)
    d = m*n;
    Tmn = falses(d,d);
    i = 1:d;
    rI = 1 .+ m.*(i.-1)-(m*n-1).*((i.-1).÷n);
    Tmn[[CartesianIndex(v,w) for (w,v) in zip(rI,i)]] .= 1;
    return Tmn
end


text:
const FLAG_SIMULATION_BERNOULLI_VARY_RANK_SMALL                 = false
const FLAG_SIMULATION_BERNOULLI_VARY_RANK_MEDIUM                = true
const FLAG_SIMULATION_BERNOULLI_VARY_RANK_LARGE                 = false

const FLAG_SIMULATION_BERNOULLI_VARY_MISSING_PERCENTAGE_SMALL   = false
const FLAG_SIMULATION_BERNOULLI_VARY_MISSING_PERCENTAGE_MEDIUM  = false
const FLAG_SIMULATION_BERNOULLI_VARY_MISSING_PERCENTAGE_LARGE   = false

FLAG_SIMULATION_BERNOULLI_VARY_RANK_SMALL ?
  include("simulation_bernoulli_vary_rank_small.jl")  : nothing

FLAG_SIMULATION_BERNOULLI_VARY_RANK_MEDIUM ?
  include("simulation_bernoulli_vary_rank_medium.jl") : nothing

FLAG_SIMULATION_BERNOULLI_VARY_RANK_LARGE ?
  include("simulation_bernoulli_vary_rank_large.jl") : nothing

# @testset "$(format("ADMM Algorithm: Small Input Simulation [Bernoulli 400 x 400]"))" begin
#   if SIMULATION_STATUS_BERNOULLI_VARY_RANK_SMALL == false
#     let
#       Random.seed!(65536)
#       ROW = 400
#       COL = 400
#       for input_rank in 1:2:400
#         @printf("small case: rank = %d\n", input_rank)
#         dist = Bernoulli()
#         timer = TimerOutput()
#         Base.Filesystem.mkpath("./test_result/bernoulli/small(400x400)")
#         io = open("./test_result/bernoulli/small(400x400)/rank"*string(input_rank)*".log", "w")
#         truth_matrix        = rand([(FixedRankMatrix(Distributions.Bernoulli(0.5), rank = input_rank), ROW, COL)])
#         sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
#         input_matrix        = sample_model.draw(truth_matrix)
#         manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
#         manual_type_matrix .= :Bernoulli
#         @timeit timer  "Bernoulli(400x400)" * "| rank="*string(input_rank) begin
#           completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
#                                                              maxiter         = 200,
#                                                              ρ               = 0.3,
#                                                              use_autodiff    = false,
#                                                              gd_iter         = 3,
#                                                              debug_mode      = false,
#                                                              project_rank    = nothing,
#                                                              io              = io,
#                                                              type_assignment = manual_type_matrix)
#         end
#         log_simulation_result(Bernoulli(), completed_matrix, truth_matrix, type_tracker,tracker, io = io)
#         show(io, MIME("text/plain"), timer)
#       end
#     end
#   end
# end

# @testset "$(format("ADMM Algorithm: Medium Input Simulation [Bernoulli 2000 x 2000]"))" begin
#   if SIMULATION_STATUS_BERNOULLI_VARY_RANK_MEDIUM == false
#     let
#       Random.seed!(65536)
#       ROW = 2000
#       COL = 2000
#       # for input_rank in union(1, collect(10:10:500))
#       for input_rank in collect(300:10:500)
#         @printf("medium case: rank = %d\n", input_rank)
#         dist = Bernoulli()
#         timer = TimerOutput()
#         RESULTS_DIR    = GLOBAL_SIMULATION_RESULTS_DIR * "bernoulli/medium(2000x2000)(vary_rank)/" * "rank" * string(input_rank) * "/"
#         LOG_FILE_NAME  = "io.log"
#         DATA_FILE_NAME = "saved_variables.h5"
#         LOG_FILE_PATH  = RESULTS_DIR * LOG_FILE_NAME
#         DATA_FILE_PATH = RESULTS_DIR * DATA_FILE_NAME
#         Base.Filesystem.mkpath(RESULTS_DIR)
#         io = open(LOG_FILE_PATH, "w")
#         truth_matrix        = rand([(FixedRankMatrix(Distributions.Bernoulli(0.5), rank = input_rank), ROW, COL)])
#         sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
#         input_matrix        = sample_model.draw(truth_matrix)
#         manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
#         manual_type_matrix .= :Bernoulli
#         @timeit timer  "Bernoulli(2000x2000)" * "| rank=" * string(input_rank) begin
#           completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
#                                                              maxiter         = 200,
#                                                               ρ               = 0.3,
#                                                              use_autodiff    = false,
#                                                              gd_iter         = 3,
#                                                              debug_mode      = false,
#                                                              project_rank    = input_rank,
#                                                              io              = io,
#                                                              type_assignment = manual_type_matrix)
#         end

#         predicted_matrix = predict(MatrixCompletionModel(),
#                                    completed_matrix = completed_matrix,
#                                    type_tracker     = type_tracker)

#         summary_object   = summary(MatrixCompletionModel(),
#                                    predicted_matrix = predicted_matrix,
#                                    truth_matrix     = truth_matrix,
#                                    type_tracker     = type_tracker,
#                                    tracker          = tracker)
#         pickle(DATA_FILE_PATH,
#                "missing_idx"      => type_tracker[:Missing],
#                "completed_matrix" => completed_matrix,
#                "predicted_matrix" => predicted_matrix,
#                "truth_matrix"     => truth_matrix,
#                "summary"          => summary_object)

#         log_simulation_result(Bernoulli(), completed_matrix, truth_matrix, type_tracker,tracker, io = io)
#         print(io, JSON.json(summary_object, 4))
#         print(io, timer)
#         close(io)
#       end
#     end
#   else
#     @info("Already Completed Simualtion[Bernoulli vary rank][MEDIUM]")
#   end
# end

# @testset "$(format("ADMM Algorithm: Medium Input Simulation [Bernoulli 2000 x 2000]"))" begin
#   if SIMULATION_STATUS_BERNOULLI_VARY_RANK_MEDIUM == false
#     let
#       Random.seed!(65536)
#       ROW = 100
#       COL = 100
#       for input_rank in union(1, collect(10:10:500))
#         @printf("medium case: rank = %d\n", input_rank)
#         dist = Bernoulli()
#         timer = TimerOutput()
#         # Base.Filesystem.mkpath("./test_result/bernoulli/medium(2000x2000)")
#         # io = open("./test_result/bernoulli/medium(2000x2000)/rank"*string(input_rank)*".log", "w")
#         io = Base.stdout
#         truth_matrix        = rand([(FixedRankMatrix(Distributions.Bernoulli(0.5), rank = input_rank), ROW, COL)])
#         sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
#         input_matrix        = sample_model.draw(truth_matrix)
#         manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
#         manual_type_matrix .= :Bernoulli
#         @timeit timer  "Bernoulli(2000x2000)" * "| rank="*string(input_rank) begin
#           completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
#                                                              maxiter         = 200,
#                                                              ρ               = 0.3,
#                                                              use_autodiff    = false,
#                                                              gd_iter         = 3,
#                                                              debug_mode      = false,
#                                                              project_rank    = input_rank,
#                                                              io              = io,
#                                                              type_assignment = manual_type_matrix)
#         end
#         log_simulation_result(Bernoulli(), completed_matrix, truth_matrix, type_tracker,tracker, io = io)
#         show(io, MIME("text/plain"), timer)
#       end
#     end
#   end
# end


text:
module chapters_deployment
using Dash, DashCoreComponents, DashHtmlComponents

include("../../utils.jl")

app = dash()

app.layout = html_div() do
    html_h1("Deploying Dash Apps"),
    dcc_markdown("
By default, Dash apps run on `localhost` - you can only access them on your
own machine. To share a Dash app, you need to 'deploy' your Dash app to a
server and open up the server's firewall to the public or to a restricted
set of IP addresses.

## Dash Enterprise

[Dash Enterprise](https://plotly.com/dash/?_ga=2.249471751.1080104966.1578062860-1986131108.1567098614) is Plotly's commercial product for deploying Dash
Apps on your company's servers or on AWS, Google Cloud, or Azure. It
offers an enterprise-wide Dash App Portal, easy git-based deployment,
automatic URL namespacing, built-in SSL support, LDAP authentication, and
more. [Learn more about Dash Enterprise](https://plotly.com/dash?_ga=2.249471751.1080104966.1578062860-1986131108.1567098614) or [get in touch to start a
trial](https://plotly.com/get-demo/?_ga=2.48144423.1080104966.1578062860-1986131108.1567098614).

For existing customers, see the [Dash Enterprise Documentation](https://dash.plotly.com/dash-enterprise).

## Dash and HTTP.jl

Dash apps are web applications. Dash relies upon HTTP.jl to provide client
and server functionality. When you create a Dash for Julia application, the
resulting object provides an interface that wraps this functionality:

```julia
    using Dash

    app = dash()
```

## Heroku Example

Heroku is one of the easiest platforms for deploying and managing public web applications.

Here is a simple example. This example requires a Heroku account and `git`.

---

**Step 1.** Create a new folder for your project:

```
    \$ mkdir dash_app_example
    \$ cd dash_app_example
```

---

**Step 2.** Initialize the folder with `git`

```
    \$ git init        # initializes an empty git repo
```

---

**Step 3.** Initialize the folder with a sample app (`app.jl`), and a `.gitignore` file (not required, but will avoid committing any files that aren't necessary for your app to function).

Add the following Julia application to your project folder:

**`app.jl`**

```julia
    using Dash, DashHtmlComponents, DashCoreComponents;

    app = dash();

    app.layout = html_div() do
        dcc_dropdown(
            id = \"dropdown\",
            options = [
                Dict(\"label\" => \"Los Angeles\", \"value\" => \"LA\"),
                Dict(\"label\" => \"New York City\", \"value\" => \"NYC\"),
                Dict(\"label\" => \"Montreal\", \"value\" => \"MTL\")
            ],
            value = \"LA\"),
        html_div(id = \"display-value\")
    end;

    callback!(app,
        Output(\"display-value\", \"children\"),
        Input(\"dropdown\", \"value\")) do value
            \"You've entered \$(value)\"
    end;

    # the default Dash port is 8050, but for Heroku deployments, it's
    # important to ensure the server port matches the PORT environment
    # variable; when running the app on your own machine, omit this
    # unless PORT is set.
    port = parse(Int64, ENV[\"PORT\"]);

    run_server(app, \"0.0.0.0\", port)
```

---

**Step 4.** Use Julia's package manager to generate a file (`Project.toml`) containing all the dependencies for your application.

Make note of any packages you will need for your application once deployed, and install them using `Pkg`. For example, given
our application above, we'll want to launch Julia and then type `]` at the prompt. You should see a new prompt, which includes
the current Julia version and `pkg>`.

From within the directory you created to store your new Dash application, type `activate .` to create a new environment for
this deployment:

```julia
    (@v1.5) pkg> activate .
    Activating new environment at `/dash-projects/dash_app_example/Project.toml`
```

Now, from the same prompt, install the packages (edit as needed for your own application):

```julia
   (dash_app_example) pkg> add Dash, DashHtmlComponents, DashCoreComponents
```

Julia should now display a message indicating that the registry is being updated, and that package versions are being resolved
and installed. Once finished, you should see the `pkg>` prompt again. Hit the backspace key (or hold down the control and C keys), and this should return
you to the `julia>` prompt. From here we can exit Julia, and we'll include the `Project.toml` and `Manifest.toml` files that
were generated in Step 5.

---

**Step 5.** Initialize Heroku, add files to Git, and deploy

```
    \$ heroku create my-dash-app --buildpack https://github.com/Optomatica/heroku-buildpack-julia # change my-dash-app to a unique name
    \$ git add . # add all files to git
    \$ git commit -m 'Initial app boilerplate'
    \$ git push heroku master # deploy code to Heroku
    \$ heroku ps:scale web=1  # run the app with one Heroku 'dyno'
```

You should be able to access your app at `https://my-dash-app.herokuapp.com` (changing my-dash-app to the name of your app).

---

**Step 6.** Update the code and redeploy

When you modify app.jl with your own code, you will need to add the changes to git and push those changes to heroku.

```
    \$ git status # view the changes
    \$ git add .  # add all the changes
    \$ git commit -m 'a description of the changes'
    \$ git push heroku master
```

If you're ready to take your apps to the next level, and deliver interactive analytics at scale, we invite you to learn more about Dash Enterprise.
[Click here for more information](https://plotly.com/dash/?_ga=2.176345125.1075922756.1562168385-916141078.1562168385) or [get in touch](https://plotly.typeform.com/to/rkO85m?_ga=2.176345125.1075922756.1562168385-916141078.1562168385).
    "),

    html_hr(),
    dcc_markdown("[Back to the table of contents](/)")
end

end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
using BufferedStream
using Base.Test

# Construct a stream and a view of said stream
stream = LinkedStream{Float64}()
view = StreamView(stream)

# Initially, any kind of access on view is going to throw
@test_throws BoundsError view[-1]
@test_throws BoundsError view[0]
@test_throws BoundsError view[1]
@test view.chunk == stream.latest
@test isnull(view.chunk.next)

# Add a chunk of data
data = Float64[3.1,4.1,5.9,2.6,5.3]
push!(stream, data)
@test view.chunk != stream.latest
@test view.chunk.next.value == stream.latest
@test isnull(stream.latest.next)

# Print then out to see what they look like:
display(stream)
println()
display(view)
println()
display(view.chunk)
println()

# Ensure that the data came through okay, and that out of bounds accesses are recognized as such
@test length(view) == length(data)
for idx in 1:length(data)
    @test view[idx] == data[idx]
end
@test_throws BoundsError view[-1]
@test_throws BoundsError view[0]
@test_throws BoundsError view[6]
@test_throws DomainError shift_origin!(view, -1)

# Verify that the LinkedStream does not copy data:
data[1] = 3.2
@test view[1] == 3.2
data[1] = 3.1
@test view[1] == 3.1

# Add a second chunk, ensure that it is seamlessly integrated into our view
push!(stream, data)
for idx in 1:length(data)
    @test view[idx] == data[idx]
    @test view[idx + length(data)] == data[idx]
end

# Move view forward a few samples, ensure the shift worked:
shift_origin!(view, 2)
@test length(view) == 2*length(data) - 2
for idx in 1:length(view)
    @test view[idx] == data[mod1(idx + 2, length(data))]
end

# Test that the underlying data is still the original representations
@test view.chunk.data == data

# Move view forward enough to get into the next chunk
shift_origin!(view, 4)
@test length(view) == length(data) - 1
for idx in 1:length(view)
    @test view[idx] == data[idx + 1]
end

# Now add a bunch of chunks
push!(stream, Float64[1:100 ...])
push!(stream, Float64[1:100 ...])
push!(stream, Float64[1:100 ...])
push!(stream, Float64[1:100 ...])

# Ensure our length can still be found, despite all appearances to the contrary
@test length(view) == 404

# Shift the view, and start playing with more complex operations/indexing
shift_origin!(view, 157)
@test view[48:147] == Float64[1:100 ...]
@test length(view[48:end]) == 200
@test mean(view[48:end]) == mean(1:100)
@test sum(view[48:57]) == sum(1:10)

# Create a second view off of this view, and a third view off of the stream
view2 = StreamView(view)
view3 = StreamView(stream)

# Test to make sure that the linked list looks like we'd expect
@test view3.sample_offset == 100
@test view2.sample_offset == view.sample_offset
@test view2.chunk == view.chunk
@test view2.chunk.next.value.next.value == view3.chunk
@test isnull(view3.chunk.next)

# Test to make sure the values look like we'd expect
@test view2.abs_idx == view.abs_idx
@test view3.abs_idx == view.abs_idx + 247
@test view2[1] == view[1]
@test_throws BoundsError view3[1]