text:
#=

https://publiclab.org/notes/cfastie/08-26-2014/new-ndvi-colormap

This new colormap has several useful features (for about 90%
of us).

ColorSchemes's values are [0..1], the following description
talks about values in the range [-1..+1].

1 The left half of the colormap has three gradients between
black and white. This does not allow you to discern the
value from the tone, but the values of NDVI below zero
indicate zero photosynthesis, so we generally do not care
exactly what the value is. The multiple gradients preserve
the detail of non-plants in the NDVI image. It just makes
the NDVI image easier to look at because you can recognize
objects and textures that are not foliage.

2 The boundary between grayscale and color is not at zero.
Live foliage generally does not have NDVI values below 0.1,
so that is the boundary between grayscale and color. With
this colormap, anything in grayscale is probably not a
plant. This allows a more precise differentiation between
plant and non-plant when the NDVI values are calibrated.

3 There is a narrow band of violet between 0.1 and 0.2 which
could represent very low photosynthetic activity, but might
also be noise or error.

4 The primary gradient of photosynthetic activity is from
NDVI values from 0.2 to 0.9, and that is represented with a
classic heat map from green to yellow to red. It's a little
bit counter intuitive because green does not represent the
healthiest plants, but the heat map metaphor seems to work
well for most people.

5 The highest values (> 0.9) are colored magenta. Foliage
generally does not have NDVI values this high, so this color
represents non-plants. I did not make it a gray so it can be
distinguished from low values. In many cases, DIY NDVI
values above 0.9 are artifacts where the image is very dark
or very bright.

=#

loadcolorscheme(:fastie, [
RGB(1.0, 1.0, 1.0),
RGB(0.9803921568627451, 0.9803921568627451, 0.9803921568627451),
RGB(0.9647058823529412, 0.9647058823529412, 0.9647058823529412),
RGB(0.9490196078431372, 0.9490196078431372, 0.9490196078431372),
RGB(0.9333333333333333, 0.9333333333333333, 0.9333333333333333),
RGB(0.9137254901960784, 0.9137254901960784, 0.9137254901960784),
RGB(0.8980392156862745, 0.8980392156862745, 0.8980392156862745),
RGB(0.8823529411764706, 0.8823529411764706, 0.8823529411764706),
RGB(0.8666666666666667, 0.8666666666666667, 0.8666666666666667),
RGB(0.8470588235294118, 0.8470588235294118, 0.8470588235294118),
RGB(0.8313725490196079, 0.8313725490196079, 0.8313725490196079),
RGB(0.8156862745098039, 0.8156862745098039, 0.8156862745098039),
RGB(0.8, 0.8, 0.8),
RGB(0.7843137254901961, 0.7843137254901961, 0.7843137254901961),
RGB(0.7647058823529411, 0.7647058823529411, 0.7647058823529411),
RGB(0.7490196078431373, 0.7490196078431373, 0.7490196078431373),
RGB(0.7333333333333333, 0.7333333333333333, 0.7333333333333333),
RGB(0.7176470588235294, 0.7176470588235294, 0.7176470588235294),
RGB(0.6980392156862745, 0.6980392156862745, 0.6980392156862745),
RGB(0.6823529411764706, 0.6823529411764706, 0.6823529411764706),
RGB(0.6666666666666666, 0.6666666666666666, 0.6666666666666666),
RGB(0.6509803921568628, 0.6509803921568628, 0.6509803921568628),
RGB(0.6313725490196078, 0.6313725490196078, 0.6313725490196078),
RGB(0.615686274509804, 0.615686274509804, 0.615686274509804),
RGB(0.6, 0.6, 0.6),
RGB(0.5843137254901961, 0.5843137254901961, 0.5843137254901961),
RGB(0.5686274509803921, 0.5686274509803921, 0.5686274509803921),
RGB(0.5490196078431373, 0.5490196078431373, 0.5490196078431373),
RGB(0.5333333333333333, 0.5333333333333333, 0.5333333333333333),
RGB(0.5176470588235295, 0.5176470588235295, 0.5176470588235295),
RGB(0.5019607843137255, 0.5019607843137255, 0.5019607843137255),
RGB(0.4823529411764706, 0.4823529411764706, 0.4823529411764706),
RGB(0.4666666666666667, 0.4666666666666667, 0.4666666666666667),
RGB(0.45098039215686275, 0.45098039215686275, 0.45098039215686275),
RGB(0.43529411764705883, 0.43529411764705883, 0.43529411764705883),
RGB(0.41568627450980394, 0.41568627450980394, 0.41568627450980394),
RGB(0.4, 0.4, 0.4),
RGB(0.3843137254901961, 0.3843137254901961, 0.3843137254901961),
RGB(0.3686274509803922, 0.3686274509803922, 0.3686274509803922),
RGB(0.35294117647058826, 0.35294117647058826, 0.35294117647058826),
RGB(0.3333333333333333, 0.3333333333333333, 0.3333333333333333),
RGB(0.3176470588235294, 0.3176470588235294, 0.3176470588235294),
RGB(0.30196078431372547, 0.30196078431372547, 0.30196078431372547),
RGB(0.28627450980392155, 0.28627450980392155, 0.28627450980392155),
RGB(0.26666666666666666, 0.26666666666666666, 0.26666666666666666),
RGB(0.25098039215686274, 0.25098039215686274, 0.25098039215686274),
RGB(0.23529411764705882, 0.23529411764705882, 0.23529411764705882),
RGB(0.2196078431372549, 0.2196078431372549, 0.2196078431372549),
RGB(0.20392156862745098, 0.20392156862745098, 0.20392156862745098),
RGB(0.2196078431372549, 0.2196078431372549, 0.2196078431372549),
RGB(0.23529411764705882, 0.23529411764705882, 0.23529411764705882),
RGB(0.25098039215686274, 0.25098039215686274, 0.25098039215686274),
RGB(0.26666666666666666, 0.26666666666666666, 0.26666666666666666),
RGB(0.28627450980392155, 0.28627450980392155, 0.28627450980392155),
RGB(0.30196078431372547, 0.30196078431372547, 0.30196078431372547),
RGB(0.3176470588235294, 0.3176470588235294, 0.3176470588235294),
RGB(0.3333333333333333, 0.3333333333333333, 0.3333333333333333),
RGB(0.35294117647058826, 0.35294117647058826, 0.35294117647058826),
RGB(0.3686274509803922, 0.3686274509803922, 0.3686274509803922),
RGB(0.3843137254901961, 0.3843137254901961, 0.3843137254901961),
RGB(0.4, 0.4, 0.4),
RGB(0.41568627450980394, 0.41568627450980394, 0.41568627450980394),
RGB(0.43529411764705883, 0.43529411764705883, 0.43529411764705883),
RGB(0.45098039215686275, 0.45098039215686275, 0.45098039215686275),
RGB(0.4666666666666667, 0.4666666666666667, 0.4666666666666667),
RGB(0.4823529411764706, 0.4823529411764706, 0.4823529411764706),
RGB(0.5019607843137255, 0.5019607843137255, 0.5019607843137255),
RGB(0.5176470588235295, 0.5176470588235295, 0.5176470588235295),
RGB(0.5333333333333333, 0.5333333333333333, 0.5333333333333333),
RGB(0.5490196078431373, 0.5490196078431373, 0.5490196078431373),
RGB(0.5686274509803921, 0.5686274509803921, 0.5686274509803921),
RGB(0.5843137254901961, 0.5843137254901961, 0.5843137254901961),
RGB(0.6, 0.6, 0.6),
RGB(0.615686274509804, 0.615686274509804, 0.615686274509804),
RGB(0.6313725490196078, 0.6313725490196078, 0.6313725490196078),
RGB(0.6509803921568628, 0.6509803921568628, 0.6509803921568628),
RGB(0.6666666666666666, 0.6666666666666666, 0.6666666666666666),
RGB(0.6823529411764706, 0.6823529411764706, 0.6823529411764706),
RGB(0.6980392156862745, 0.6980392156862745, 0.6980392156862745),
RGB(0.7176470588235294, 0.7176470588235294, 0.7176470588235294),
RGB(0.7333333333333333, 0.7333333333333333, 0.7333333333333333),
RGB(0.7490196078431373, 0.7490196078431373, 0.7490196078431373),
RGB(0.7647058823529411, 0.7647058823529411, 0.7647058823529411),
RGB(0.7843137254901961, 0.7843137254901961, 0.7843137254901961),
RGB(0.8, 0.8, 0.8),
RGB(0.8156862745098039, 0.8156862745098039, 0.8156862745098039),
RGB(0.8313725490196079, 0.8313725490196079, 0.8313725490196079),
RGB(0.8470588235294118, 0.8470588235294118, 0.8470588235294118),
RGB(0.8666666666666667, 0.8666666666666667, 0.8666666666666667),
RGB(0.8823529411764706, 0.8823529411764706, 0.8823529411764706),
RGB(0.8980392156862745, 0.8980392156862745, 0.8980392156862745),
RGB(0.9137254901960784, 0.9137254901960784, 0.9137254901960784),
RGB(0.9333333333333333, 0.9333333333333333, 0.9333333333333333),
RGB(0.9490196078431372, 0.9490196078431372, 0.9490196078431372),
RGB(0.9647058823529412, 0.9647058823529412, 0.9647058823529412),
RGB(0.9803921568627451, 0.9803921568627451, 0.9803921568627451),
RGB(1.0, 1.0, 1.0),
RGB(0.9803921568627451, 0.9803921568627451, 0.9803921568627451),
RGB(0.9607843137254902, 0.9607843137254902, 0.9607843137254902),
RGB(0.9411764705882353, 0.9411764705882353, 0.9411764705882353),
RGB(0.9215686274509803, 0.9215686274509803, 0.9215686274509803),
RGB(0.9019607843137255, 0.9019607843137255, 0.9019607843137255),
RGB(0.8823529411764706, 0.8823529411764706, 0.8823529411764706),
RGB(0.8627450980392157, 0.8627450980392157, 0.8627450980392157),
RGB(0.8431372549019608, 0.8431372549019608, 0.8431372549019608),
RGB(0.8235294117647058, 0.8235294117647058, 0.8235294117647058),
RGB(0.803921568627451, 0.803921568627451, 0.803921568627451),
RGB(0.7843137254901961, 0.7843137254901961, 0.7843137254901961),
RGB(0.7647058823529411, 0.7647058823529411, 0.7647058823529411),
RGB(0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
RGB(0.7254901960784313, 0.7254901960784313, 0.7254901960784313),
RGB(0.7058823529411765, 0.7058823529411765, 0.7058823529411765),
RGB(0.6862745098039216, 0.6862745098039216, 0.6862745098039216),
RGB(0.6666666666666666, 0.6666666666666666, 0.6666666666666666),
RGB(0.6470588235294118, 0.6470588235294118, 0.6470588235294118),
RGB(0.6274509803921569, 0.6274509803921569, 0.6274509803921569),
RGB(0.6078431372549019, 0.6078431372549019, 0.6078431372549019),
RGB(0.592156862745098, 0.592156862745098, 0.592156862745098),
RGB(0.5725490196078431, 0.5725490196078431, 0.5725490196078431),
RGB(0.5529411764705883, 0.5529411764705883, 0.5529411764705883),
RGB(0.5333333333333333, 0.5333333333333333, 0.5333333333333333),
RGB(0.5137254901960784, 0.5137254901960784, 0.5137254901960784),
RGB(0.49411764705882355, 0.49411764705882355, 0.49411764705882355),
RGB(0.4745098039215686, 0.4745098039215686, 0.4745098039215686),
RGB(0.4549019607843137, 0.4549019607843137, 0.4549019607843137),
RGB(0.43529411764705883, 0.43529411764705883, 0.43529411764705883),
RGB(0.41568627450980394, 0.41568627450980394, 0.41568627450980394),
RGB(0.396078431372549, 0.396078431372549, 0.396078431372549),
RGB(0.3764705882352941, 0.3764705882352941, 0.3764705882352941),
RGB(0.3568627450980392, 0.3568627450980392, 0.3568627450980392),
RGB(0.33725490196078434, 0.33725490196078434, 0.33725490196078434),
RGB(0.3176470588235294, 0.3176470588235294, 0.3176470588235294),
RGB(0.2980392156862745, 0.2980392156862745, 0.2980392156862745),
RGB(0.2784313725490196, 0.2784313725490196, 0.2784313725490196),
RGB(0.25882352941176473, 0.25882352941176473, 0.25882352941176473),
RGB(0.23921568627450981, 0.23921568627450981, 0.23921568627450981),
RGB(0.2196078431372549, 0.2196078431372549, 0.2196078431372549),
RGB(0.25882352941176473, 0.25882352941176473, 0.3137254901960784),
RGB(0.30196078431372547, 0.30196078431372547, 0.4117647058823529),
RGB(0.3411764705882353, 0.3411764705882353, 0.5098039215686274),
RGB(0.3843137254901961, 0.3843137254901961, 0.6078431372549019),
RGB(0.4235294117647059, 0.4235294117647059, 0.7058823529411765),
RGB(0.4666666666666667, 0.4666666666666667, 0.803921568627451),
RGB(0.5058823529411764, 0.5058823529411764, 0.9019607843137255),
RGB(0.5490196078431373, 0.5490196078431373, 1.0),
RGB(0.5137254901960784, 0.5764705882352941, 0.9372549019607843),
RGB(0.47843137254901963, 0.6039215686274509, 0.8745098039215686),
RGB(0.44313725490196076, 0.6313725490196078, 0.8117647058823529),
RGB(0.4117647058823529, 0.6588235294117647, 0.7490196078431373),
RGB(0.3764705882352941, 0.6862745098039216, 0.6862745098039216),
RGB(0.3411764705882353, 0.7176470588235294, 0.6235294117647059),
RGB(0.3058823529411765, 0.7450980392156863, 0.5607843137254902),
RGB(0.27450980392156865, 0.7725490196078432, 0.4980392156862745),
RGB(0.23921568627450981, 0.8, 0.43529411764705883),
RGB(0.20392156862745098, 0.8274509803921568, 0.37254901960784315),
RGB(0.16862745098039217, 0.8588235294117647, 0.30980392156862746),
RGB(0.13725490196078433, 0.8862745098039215, 0.24705882352941178),
RGB(0.10196078431372549, 0.9137254901960784, 0.1843137254901961),
RGB(0.06666666666666667, 0.9411764705882353, 0.12156862745098039),
RGB(0.03137254901960784, 0.9686274509803922, 0.058823529411764705),
RGB(0.0, 1.0, 0.0),
RGB(0.027450980392156862, 1.0, 0.0),
RGB(0.058823529411764705, 1.0, 0.0),
RGB(0.09019607843137255, 1.0, 0.0),
RGB(0.12156862745098039, 1.0, 0.0),
RGB(0.15294117647058825, 1.0, 0.0),
RGB(0.1843137254901961, 1.0, 0.0),
RGB(0.21568627450980393, 1.0, 0.0),
RGB(0.24705882352941178, 1.0, 0.0),
RGB(0.2784313725490196, 1.0, 0.0),
RGB(0.30980392156862746, 1.0, 0.0),
RGB(0.3411764705882353, 1.0, 0.0),
RGB(0.37254901960784315, 1.0, 0.0),
RGB(0.403921568627451, 1.0, 0.0),
RGB(0.43529411764705883, 1.0, 0.0),
RGB(0.4666666666666667, 1.0, 0.0),
RGB(0.4980392156862745, 1.0, 0.0),
RGB(0.5294117647058824, 1.0, 0.0),
RGB(0.5607843137254902, 1.0, 0.0),
RGB(0.592156862745098, 1.0, 0.0),
RGB(0.6235294117647059, 1.0, 0.0),
RGB(0.6549019607843137, 1.0, 0.0),
RGB(0.6862745098039216, 1.0, 0.0),
RGB(0.7176470588235294, 1.0, 0.0),
RGB(0.7490196078431373, 1.0, 0.0),
RGB(0.7803921568627451, 1.0, 0.0),
RGB(0.8117647058823529, 1.0, 0.0),
RGB(0.8431372549019608, 1.0, 0.0),
RGB(0.8745098039215686, 1.0, 0.0),
RGB(0.9058823529411765, 1.0, 0.0),
RGB(0.9372549019607843, 1.0, 0.0),
RGB(0.9686274509803922, 1.0, 0.0),
RGB(1.0, 1.0, 0.0),
RGB(1.0, 0.9764705882352941, 0.0),
RGB(1.0, 0.9568627450980393, 0.0),
RGB(1.0, 0.9372549019607843, 0.0),
RGB(1.0, 0.9137254901960784, 0.0),
RGB(1.0, 0.8941176470588236, 0.0),
RGB(1.0, 0.8745098039215686, 0.0),
RGB(1.0, 0.8509803921568627, 0.0),
RGB(1.0, 0.8313725490196079, 0.0),
RGB(1.0, 0.8117647058823529, 0.0),
RGB(1.0, 0.788235294117647, 0.0),
RGB(1.0, 0.7686274509803922, 0.0),
RGB(1.0, 0.7490196078431373, 0.0),
RGB(1.0, 0.7254901960784313, 0.0),
RGB(1.0, 0.7058823529411765, 0.0),
RGB(1.0, 0.6862745098039216, 0.0),
RGB(1.0, 0.6666666666666666, 0.0),
RGB(1.0, 0.6431372549019608, 0.0),
RGB(1.0, 0.6235294117647059, 0.0),
RGB(1.0, 0.6039215686274509, 0.0),
RGB(1.0, 0.5803921568627451, 0.0),
RGB(1.0, 0.5607843137254902, 0.0),
RGB(1.0, 0.5411764705882353, 0.0),
RGB(1.0, 0.5176470588235295, 0.0),
RGB(1.0, 0.4980392156862745, 0.0),
RGB(1.0, 0.47843137254901963, 0.0),
RGB(1.0, 0.4549019607843137, 0.0),
RGB(1.0, 0.43529411764705883, 0.0),
RGB(1.0, 0.41568627450980394, 0.0),
RGB(1.0, 0.39215686274509803, 0.0),
RGB(1.0, 0.37254901960784315, 0.0),
RGB(1.0, 0.35294117647058826, 0.0),
RGB(1.0, 0.3333333333333333, 0.0),
RGB(1.0, 0.30980392156862746, 0.0),
RGB(1.0, 0.2901960784313726, 0.0),
RGB(1.0, 0.27058823529411763, 0.0),
RGB(1.0, 0.24705882352941178, 0.0),
RGB(1.0, 0.22745098039215686, 0.0),
RGB(1.0, 0.20784313725490197, 0.0),
RGB(1.0, 0.1843137254901961, 0.0),
RGB(1.0, 0.16470588235294117, 0.0),
RGB(1.0, 0.1450980392156863, 0.0),
RGB(1.0, 0.12156862745098039, 0.0),
RGB(1.0, 0.10196078431372549, 0.0),
RGB(1.0, 0.08235294117647059, 0.0),
RGB(1.0, 0.058823529411764705, 0.0),
RGB(1.0, 0.0392156862745098, 0.0),
RGB(1.0, 0.0196078431372549, 0.0),
RGB(1.0, 0.0, 0.0),
RGB(1.0, 0.0, 0.058823529411764705),
RGB(1.0, 0.0, 0.12156862745098039),
RGB(1.0, 0.0, 0.1843137254901961),
RGB(1.0, 0.0, 0.24705882352941178),
RGB(1.0, 0.0, 0.30980392156862746),
RGB(1.0, 0.0, 0.37254901960784315),
RGB(1.0, 0.0, 0.43529411764705883),
RGB(1.0, 0.0, 0.4980392156862745),
RGB(1.0, 0.0, 0.5607843137254902),
RGB(1.0, 0.0, 0.6235294117647059),
RGB(1.0, 0.0, 0.6862745098039216),
RGB(1.0, 0.0, 0.7490196078431373),
RGB(1.0, 0.0, 0.8117647058823529),
RGB(1.0, 0.0, 0.8745098039215686),
RGB(1.0, 0.0, 0.9372549019607843),
], "general", "fastie, vegetation index, NDVI")

max_stars_repo_name: cormullion/ColorSchemes

max_issues_repo_name: cormullion/ColorSchemes

max_forks_repo_name: cormullion/ColorSchemes


text:
# This file is a part of BATBase.jl, licensed under the MIT License (MIT).

include("abstract_density_sample.jl")
include("abstract_density.jl")

ext: jl


text:
#
# Data structures and types for orbit maneuvers
#

"""
An abstract type for all orbit maneuvers.
"""
abstract type AbstractManeuver end

"""
An abstract type for all twobody maneuvers.
"""
abstract type TwoBodyManeuver <: AbstractManeuver end

"""
A type for constant, continuous thrust.
"""
struct ConstantManeuver <: TwoBodyManeuver

    aₜ::Unitful.Acceleration

end


text:
using Thebes, Luxor

function makemobius()
    x(u, v) = (1 + (v/2 * cos(u/2))) * cos(u)
    y(u, v) = (1 + (v/2 * cos(u/2))) * sin(u)
    z(u, v) = v/2 * sin(u/2)
    # v goes across
    # u goes round
    w = .5
    st = 2π/200
    Δ = .05
    result = []
    for u in 0:st:2π-st
        for v in -w:Δ:w
            p1 = Point3D(x(u,      v + Δ),   y(u,      v + Δ),     z(u,       v + Δ))
            p2 = Point3D(x(u + st, v + Δ),   y(u + st, v + Δ),     z(u  + st, v + Δ))
            p3 = Point3D(x(u + st, v),       y(u + st, v),         z(u + st,  v))
            p4 = Point3D(x(u,      v),       y(u,      v),         z(u,       v))
            push!(result, [p1, p2, p3, p4])
        end
    end
    return result
end

@svg begin
    background("black")
    setline(0.2)
    eyepoint(300, 300, 300)
    perspective(1500)
    mb = makemobius()
    for p in mb
        sethue([Luxor.julia_green, Luxor.julia_red,Luxor.julia_purple, Luxor.julia_blue][rand(1:end)])
        pin(100p, gfunction  = (p3l, p2l) -> begin
            setopacity(0.85)
            sethue("grey20")
            poly(p2l, :fillpreserve, close=true)
            setopacity(1)
            sethue("white")
            strokepath()
        end)
    end
end 1200 1200 "/tmp/mo.svg"


text:
function overlap( p, q, ϵ )

   for r in q
       all(abs.(p .- r) .< 2ϵ) && (return true)
   end

   return false

end

function box_particles(rng, n, ϵ)

    l1 = 0.5
    l2 = 0.5

    cost = 1 / sqrt(2)
    J4 = [cost cost; -cost cost]

    q = Vector{Float64}[]
    push!(q, [ 0,           0])
    push!(q, [ 0.0685359 ,  0.15134794])
    push!(q, [-0.37969185, -0.02715423])
    push!(q, [-0.23362135,  0.07949096])
    push!(q, [-0.03214348,  0.07088258])
    push!(q, [ 0.04120227,  0.33615893])
    push!(q, [-0.1507703 ,  0.32089688])
    push!(q, [-0.2       , -0.2])
    push!(q, [ 0.25      , -0.05])
    push!(q, [ 0.1       , -0.2])

    p = zeros(2)

    for klm in 11:n
        frein = 0
        overlap = 1
        while (overlap == 1) && frein < 10000

            p0 = [2 * l1 * rand(rng) - l1, 2 * l2 * rand(rng) - l2]
            p = J4 * p0
            p = p * (1 - 4 * ϵ) / sqrt(2)
            overlap2 = 0
            for subh in 1:(klm-1)
                if norm(p - q[subh], 1) < 2ϵ
                    overlap2 = 1
                end
            end

            overlap = overlap2
            frein += 1
        end

        if frein == 10000
            print("Echec tirage initial")
        end
        push!(q, p)
    end

    vitesses = [[1, 0], [0, 1], [-1, 0], [0, -1]]
    v = [vitesses[rand(rng,1:end)] for i in 1:n]

    return q, v

end

export SquareParticles

struct SquareParticles <: Particles

    n
    q
    v
    ϵ

    function SquareParticles(q :: Vector{Vector{Float64}}, v, ϵ :: AbstractFloat)

        n = length(q)

        @assert n == length(v)

        new(n, q, v, ϵ)

    end

    function SquareParticles(rng, n, ϵ; option = :none )

        if option == :box

           q, v = box_particles(rng, n, ϵ)

        else

           q = Vector{Float64}[]

           push!(q, [0.5,0.5])

           for i in 1:n-1

               p = [0.5, 0.5]
               k = 0
               while overlap(p, q, ϵ) && k < 1000
                   rand!(rng, p)
                   k += 1
               end

               if k < 1000
                  push!(q, p)
               else
                  @error "echec du tirage $i"
               end

           end

           vitesses = [[1, 0], [0, 1], [-1, 0], [0, -1]]
           v = [vitesses[rand(rng, 1:end)] for i in 1:n]
        end

        new(n, q, v, ϵ)

    end

end

function compute_dt(p :: SquareParticles, i, j, k = 1)

    ϵ = p.ϵ
    deltat1 = Inf
    deltat2 = Inf
    x = p.q[i] .+ offset[k]
    y = p.q[j]
    v = p.v[i]
    w = p.v[j]

    if (y - x)'v > 0.0 > (y - x)'w

        dx = y[1] - x[1]
        dy = y[2] - x[2]
        adx = abs(dx)
        ady = abs(dy)
        adxy = abs(adx - ady)

        if w[1] == 1 && v[1] == -1 && ady < 2ϵ
            deltat1 = (2ϵ - dx - ady) / 2
            deltat2 = (- 2ϵ - dx + ady) / 2
        elseif v[1] == 1 && w[1] == -1 && ady < 2ϵ
            deltat1 = (2ϵ + dx - ady) / 2
            deltat2 = (- 2ϵ + dx + ady) / 2
        elseif w[2] == 1 && v[2] == -1 && adx < 2ϵ
            deltat1 = (2ϵ - adx - dy) / 2
            deltat2 = (- 2ϵ + adx - dy) / 2
        elseif v[2] == 1 && w[2] == -1 && adx < 2ϵ
            deltat1 = (2ϵ - adx + dy) / 2
            deltat2 = (- 2ϵ + adx + dy) / 2
        elseif w[1] == 1 && v[2] == 1 && adxy < 2ϵ
            deltat1 = (dy - dx - 2ϵ) / 2
            deltat2 = (dy - dx + 2ϵ) / 2
        elseif v[1] == 1 && w[2] == 1 && adxy < 2ϵ
            deltat1 = (-dy + dx - 2ϵ) / 2
            deltat2 = (-dy + dx + 2ϵ) / 2
        elseif w[1] == 1 && v[2] == -1 && adxy < 2ϵ
            deltat1 = (- dx - dy - 2ϵ) / 2
            deltat2 = (- dx - dy + 2ϵ) / 2
        elseif v[1] == 1 && w[2] == -1 && adxy < 2ϵ
            deltat1 = (dx + dy - 2ϵ) / 2
            deltat2 = (dx + dy + 2ϵ) / 2
        elseif w[1] == -1 && v[2] == 1 && adxy < 2ϵ
            deltat1 = (dx + dy + 2ϵ) / 2
            deltat2 = (dx + dy - 2ϵ) / 2
        elseif v[1] == -1 && w[2] == 1 && adxy < 2ϵ
            deltat1 = (- dx - dy + 2ϵ) / 2
            deltat2 = (- dx - dy - 2ϵ) / 2
        elseif w[1] == -1 && v[2] == -1 && adxy < 2ϵ
            deltat1 = (dx - dy + 2ϵ) / 2
            deltat2 = (dx - dy - 2ϵ) / 2
        elseif v[1] == -1 && w[2] == -1 && adxy < 2ϵ
            deltat1 = (- dx + dy + 2ϵ) / 2
            deltat2 = (- dx + dy - 2ϵ) / 2
        end

    end

    return min(deltat1, deltat2)

end

function compute_dt(p :: SquareParticles, j)

    ϵ = p.ϵ
    t = Inf
    x = p.q[j]
    v = p.v[j]
    i = 1

    if x[2] > 0 && v[1] == 1
        t = (0.5 - (x[1] + x[2]) - ϵ) / (v[1] + v[2])
        i = 1
    elseif x[1] > 0 && v[2] == 1
        t = (0.5 - (x[1] + x[2]) - ϵ) / (v[1] + v[2])
        i = 1
    elseif x[2] > 0 && v[1] == -1
        t = (0.5 - ϵ + (x[1] - x[2])) / (v[2] - v[1])
        i = 2
    elseif x[1] < 0 && v[2] == 1
        t = (0.5 - ϵ + (x[1] - x[2])) / (v[2] - v[1])
        i = 2
    elseif x[2] < 0 && v[1] == -1
        t = (0.5 - ϵ + x[1] + x[2]) / -(v[1] + v[2])
        i = 3
    elseif x[1] < 0 && v[2] == -1
        t = (0.5 - ϵ + x[1] + x[2]) / -(v[1] + v[2])
        i = 3
    elseif x[2] < 0 && v[1] == 1
        t = (0.5 - ϵ - (x[1] - x[2])) / (v[1] - v[2])
        i = 4
    elseif x[1] > 0 && v[2] == -1
        t = (0.5 - ϵ - (x[1] - x[2])) / (v[1] - v[2])
        i = 4
    end

    return t, i

end


text:
@testset "Restart" begin

  using Test, LinearAlgebra
  using StochasticDiffEq, DiffEqNoiseProcess
  using Random

  u₀ = [0.5]
  tstart = 0.0
  tend = 1.0
  dt = 0.005
  trange = (tstart, tend)

  f!(du,u,p,t) = du.=1.01*u
  σ!(du,u,p,t) = du.=0.87*u

  dt1 = tend/1e3
  seed = 100
  Random.seed!(seed)
  prob = SDEProblem(f!,σ!,u₀,trange)
  sol = solve(prob,EulerHeun(),dt=dt1,adaptive=false, save_noise=true, saveat=collect(trange[1]:dt1:trange[2]))

  # choose a random interval
  interval = (0.35, 0.87)

  idx1 = searchsortedfirst(sol.t, interval[1])
  idx2 = searchsortedfirst(sol.t, interval[2])

  forwardnoise = DiffEqNoiseProcess.NoiseGrid(sol.t[idx1:idx2], sol.W.W[idx1:idx2])
  forwardnoise2 = DiffEqNoiseProcess.NoiseWrapper(sol.W, indx=idx1)

  cpsol = solve(remake(prob, tspan=interval, u0=sol(interval[1]), noise=forwardnoise), sol.alg, save_noise=false; dt=dt1)
  cpsol2 = solve(remake(prob, tspan=interval, u0=sol(interval[1]), noise=forwardnoise2), sol.alg, save_noise=false; dt=dt1)

  sola = vcat(sol.u[idx1:idx2] ...)
  cpsola = vcat(cpsol.u ...)
  cpsol2a = vcat(cpsol2.u ...)

  @test isapprox(sola, cpsola, atol=1e-10 )
  @test isapprox(sola, cpsol2a, atol=1e-10 )

end


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

## dummy stub for https://github.com/JuliaBinaryWrappers/dSFMT_jll.jl

baremodule dSFMT_jll
using Base, Libdl
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libdSFMT

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libdSFMT_handle = C_NULL
libdSFMT_path = ""

if Sys.iswindows()
    const libdSFMT = "libdSFMT.dll"
elseif Sys.isapple()
    const libdSFMT = "@rpath/libdSFMT.dylib"
else
    const libdSFMT = "libdSFMT.so"
end

function __init__()
    global libdSFMT_handle = dlopen(libdSFMT)
    global libdSFMT_path = dlpath(libdSFMT_handle)
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = dirname(libdSFMT_path)
    push!(LIBPATH_list, LIBPATH[])
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libdSFMT_path() = libdSFMT_path

end  # module dSFMT_jll

max_stars_repo_name: ranocha/julia

max_issues_repo_name: ranocha/julia

max_forks_repo_name: ranocha/julia


text:
module GaussLegendre
	using LinearAlgebra
	using Printf

	function gausslegendre(n)
		# COMPUTE GAUSS-LEGENDRE NODES AND WEIGHTS USING ASYMPTOTIC EXPANSIONS.
		# COMPLEXITY O(n).

		# Nodes and weights:
		m = (n + 1) >> 1
		a = besselZeroRoots(m)
		rmul!(a, BigFloat(1) / (BigFloat(n) + BigFloat(1) / BigFloat(2)))
		x = legpts_nodes(n, a)
		w = legpts_weights(n, m, a)
		# Use symmetry to get the others:
		resize!(x, n)
		resize!(w, n)
		@inbounds for i in 1:m
			x[n + 1 - i] = x[i]
			w[n + 1 - i] = w[i]
		end
		@inbounds for i in 1:m
			x[i] = -x[i]
		end
		@inbounds mod(n, 2) == 1 && (x[m] = 0.0)
		x, w
	end

	function legpts_nodes(n, a)
		# ASYMPTOTIC EXPANSION FOR THE GAUSS-LEGENDRE NODES.
		vn = BigFloat(1) / (BigFloat(n) + BigFloat(1) / BigFloat(2))
		m = length(a)
		nodes = cot.(a)
		vn² = vn * vn
		vn⁴ = vn² * vn²
		@inbounds if n <= 255
			vn⁶ = vn⁴ * vn²
			for i in 1:m
				u = nodes[i]
				u² = u^2
				ai = a[i]
				ai² = ai * ai
				ai³ = ai² * ai
				ai⁵ = ai² * ai³
				node = ai + (u - BigFloat(1) / ai) / BigFloat(8) * vn²
				v1 = (BigFloat(6) * (BigFloat(1) + u²) / ai + BigFloat(25) / ai³ - u * muladd(BigFloat(31), u², BigFloat(33))) / BigFloat(384)
				v2 = u * @evalpoly(u², BigFloat(2595) / BigFloat(15360), BigFloat(6350) / BigFloat(15360), BigFloat(3779) / BigFloat(15360))
				v3 = (1 + u²) * (-muladd(BigFloat(31) / BigFloat(1024), u², BigFloat(11) / BigFloat(1024)) / ai +
								 u / BigFloat(512) / ai² + BigFloat(-25) / BigFloat(3072) / ai³)
				v4 = (v2 - BigFloat(1073) / BigFloat(5120) / ai⁵ + v3)
				node = muladd(v1, vn⁴, node)
				node = muladd(v4, vn⁶, node)
				nodes[i] = node
			end
		elseif n <= 3950
			for i in 1:m
				u = nodes[i]
				u² = u^2
				ai = a[i]
				ai² = ai * ai
				ai³ = ai² * ai
				node = ai + (u - 1 / ai) / 8 * vn²
				v1 = (BigFloat(6) * (BigFloat(1) + u²) / ai + BigFloat(25) / ai³ - u * muladd(BigFloat(31), u², BigFloat(33))) / BigFloat(384)
				node = muladd(v1, vn⁴, node)
				nodes[i] = node
			end
		else
			for i in 1:m
				u = nodes[i]
				ai = a[i]
				node = ai + (u - BigFloat(1) / ai) / BigFloat(8) * vn²
				nodes[i] = node
			end
		end
		@inbounds for jj = 1:m
			nodes[jj] = cos(nodes[jj])
		end
		nodes
	end

	function legpts_weights(n, m, a)
		# ASYMPTOTIC EXPANSION FOR THE GAUSS-LEGENDRE WEIGHTS.
		vn = BigFloat(1) / (BigFloat(n) + BigFloat(1) / BigFloat(2))
		vn² = vn^2
		weights = Array{Float64}(undef, m)
		if n <= 850000
			@inbounds for i in 1:m
				weights[i] = cot(a[i])
			end
		end
		# Split out the part that can be vectorized by llvm
		@inbounds if n <= 170
			for i in 1:m
				u = weights[i]
				u² = u^2
				ai = a[i]
				ai⁻¹ = 1 / ai
				ai² = ai^2
				ai⁻² = 1 / ai²
				ua = u * ai
				W1 = muladd(ua - BigFloat(1), ai⁻², BigFloat(1)) / BigFloat(8)
				W2 = @evalpoly(ai⁻², @evalpoly(u², BigFloat(-27), BigFloat(-84), BigFloat(-56)),
							   muladd(BigFloat(-3), muladd(u², BigFloat(-2), BigFloat(1)), BigFloat(6) * ua),
							   muladd(ua, BigFloat(-31), BigFloat(81))) / BigFloat(384)
				W3 = @evalpoly(ai⁻¹, @evalpoly(u², BigFloat(153) / BigFloat(1024), BigFloat(295) / BigFloat(256), BigFloat(187) / BigFloat(96),
											   BigFloat(151) / BigFloat(160)),
							   @evalpoly(u², BigFloat(-65) / BigFloat(1024), BigFloat(-119) / BigFloat(768), BigFloat(-35) / BigFloat(384)) * u,
							   @evalpoly(u², BigFloat(5) / BigFloat(512), BigFloat(15) / BigFloat(512), BigFloat(7) / BigFloat(384)),
							   muladd(u², BigFloat(1) / BigFloat(512), BigFloat(-13) / BigFloat(1536)) * u,
							   muladd(u², BigFloat(-7) / BigFloat(384), + BigFloat(53) / BigFloat(3072)),
							   BigFloat(3749) / BigFloat(15360) * u, BigFloat(-1125) / BigFloat(1024))
				weights[i] = @evalpoly(vn², BigFloat(1) / vn² + W1, W2, W3)
			end
		elseif n <= 1500
			for i in 1:m
				u = weights[i]
				u² = u^2
				ai = a[i]
				ai² = ai^2
				ai⁻² = 1 / ai²
				ua = u * ai
				W1 = muladd(ua - BigFloat(1), ai⁻², BigFloat(1)) / BigFloat(8)
				W2 = @evalpoly(ai⁻², @evalpoly(u², BigFloat(-27), BigFloat(-84), BigFloat(-56)),
							   muladd(BigFloat(-3), muladd(u², BigFloat(-2), BigFloat(1)), BigFloat(6) * ua),
							   muladd(ua, BigFloat(-31), BigFloat(81))) / BigFloat(384)
				weights[i] = muladd(vn², W2, 1 / vn² + W1)
			end
		elseif n <= 850000
			for i in 1:m
				u = weights[i]
				u² = u^2
				ai = a[i]
				ai² = ai^2
				ai⁻² = BigFloat(1) / ai²
				ua = u * ai
				W1 = muladd(ua - BigFloat(1), ai⁻², BigFloat(1)) / BigFloat(8)
				weights[i] = BigFloat(1) / vn² + W1
			end
		else
			for i in 1:m
				weights[i] = BigFloat(1) / vn²
			end
		end
		bJ1 = besselJ1(m)
		@inbounds for i in 1:m
			weights[i] = BigFloat(2) / (bJ1[i] * (a[i] / sin(a[i])) * weights[i])
		end
		return weights
	end

	const besselZeros_20 = [parse(BigFloat, "2.40482555769577244220158718235325068235397338867188"),
							parse(BigFloat, "5.52007811028631056871063265134580433368682861328125"),
							parse(BigFloat, "8.65372791291100718069628783268854022026062011718750"),
							parse(BigFloat, "11.7915344390142813324473536340519785881042480468750"),
							parse(BigFloat, "14.9309177084877866548140445956960320472717285156250"),
							parse(BigFloat, "18.0710639679109235089526919182389974594116210937500"),
							parse(BigFloat, "21.2116366298792584643706504721194505691528320312500"),
							parse(BigFloat, "24.3524715307493018201512313680723309516906738281250"),
							parse(BigFloat, "27.4934791320402531766831089043989777565002441406250"),
							parse(BigFloat, "30.6346064684319756565855641383677721023559570312500"),
							parse(BigFloat, "33.7758202135735672300143050961196422576904296875000"),
                        	parse(BigFloat, "36.9170983536640449074184289202094078063964843750000"),
							parse(BigFloat, "40.0584257646282395626258221454918384552001953125000"),
							parse(BigFloat, "43.1997917131767295018107688520103693008422851562500"),
							parse(BigFloat, "46.3411883716618149264832027256488800048828125000000"),
							parse(BigFloat, "49.4826098973978147910202096682041883468627929687500"),
							parse(BigFloat, "52.6240518411149977850982395466417074203491210937500"),
							parse(BigFloat, "55.7655107550199815591440710704773664474487304687500"),
							parse(BigFloat, "58.9069839260809402503582532517611980438232421875000"),
							parse(BigFloat, "62.0484691902267044838481524493545293807983398437500")]

	function besselZeroRoots(m)
		# BESSEL0ROOTS ROOTS OF BESSELJ(0,x). USE ASYMPTOTICS.
		# Use McMahon's expansion for the remainder (NIST, 10.21.19):
		jk = Array{BigFloat}(undef, m)
		p = (BigFloat(1071187749376) / BigFloat(315), BigFloat(0), BigFloat(-401743168) / BigFloat(105), BigFloat(0), BigFloat(120928) / BigFloat(15),
			 BigFloat(0), BigFloat(-124) / BigFloat(3), BigFloat(0), BigFloat(1), BigFloat(0))
		# First 20 are precomputed:
		@inbounds for jj = 1:min(m, 20)
			jk[jj] = besselZeros_20[jj]
		end
		@inbounds for jj = 21:min(m, 47)
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			ak82 = (BigFloat(1) / BigFloat(8) / ak)^2
			jk[jj] = ak + BigFloat(1) / BigFloat(8) / ak * @evalpoly(ak82, BigFloat(1), p[7], p[5], p[3])
		end
		@inbounds for jj = 48:min(m, 344)
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			ak82 = (BigFloat(1) / BigFloat(8) / ak)^2
			jk[jj] = ak + BigFloat(1) / BigFloat(8) / ak * @evalpoly(ak82, BigFloat(1), p[7], p[5])
		end
		@inbounds for jj = 345:min(m,13191)
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			ak82 = (BigFloat(1) / BigFloat(8) / ak)^2
			jk[jj] = ak + BigFloat(1) / BigFloat(8) / ak * muladd(ak82, p[7], BigFloat(1))
		end
		@inbounds for jj = 13192:m
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			jk[jj] = ak + BigFloat(1) / BigFloat(8) / ak
		end
		return jk
	end

	const besselJ1_10 = [parse(BigFloat, "0.26951412394191699930425603831096479552416724046623932"),
						 parse(BigFloat, "0.11578013858220369920087677406869183804800006055233685"),
						 parse(BigFloat, "0.07368635113640830090795457351166741709059792271399662"),
						 parse(BigFloat, "0.05403757319811628461998670001810067913900791658910733"),
						 parse(BigFloat, "0.04266142901724308722505899310539410056259036644129061"),
						 parse(BigFloat, "0.03524210349099609759173136127110351041151336643781974"),
						 parse(BigFloat, "0.03002107010305467407541978396539677780280158649895994"),
                     	 parse(BigFloat, "0.02614739149530809055943469810813860988716164552883349"),
						 parse(BigFloat, "0.02315912182469139499312285631116841466379136562175513"),
						 parse(BigFloat, "0.02078382912226785687257060668341061486452086057459184")]

	function besselJ1(m)
		# BESSELJ1 EVALUATE BESSELJ(1,x)^2 AT ROOTS OF BESSELJ(0,x).
		# USE ASYMPTOTICS. Use Taylor series of (NIST, 10.17.3) and McMahon's
		# expansion (NIST, 10.21.19):
		Jk2 = Array{Float64}(undef, m)
		c = (BigFloat(-171497088497) / BigFloat(15206400), BigFloat(461797) / BigFloat(1152), BigFloat(-172913) / BigFloat(8064),
			 BigFloat(151) / BigFloat(80), BigFloat(-7) / BigFloat(24), BigFloat(0), BigFloat(2))
		# First 10 are precomputed:
		@inbounds for jj = 1:min(m, 10)
			Jk2[jj] = besselJ1_10[jj]
		end
		@inbounds for jj = 11:min(m, 15)
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			ak2 = (1 / ak)^2
			Jk2[jj] = 1 / (BigFloat(π) * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3],
													  c[2], c[1]), ak2^2, c[7])
		end
		@inbounds for jj = 16:min(m, 21)
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			ak2 = (1 / ak)^2
			Jk2[jj] = 1 / (BigFloat(π) * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3], c[2]),
											ak2^2, c[7])
		end
		@inbounds for jj = 22:min(m,55)
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			ak2 = (1 / ak)^2
			Jk2[jj] = 1 / (BigFloat(π) * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3]),
											ak2^2, c[7])
		end
		@inbounds for jj = 56:min(m,279)
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			ak2 = (BigFloat(1) / ak)^2
			Jk2[jj] = BigFloat(1) / (BigFloat(π) * ak) * muladd(muladd(ak2, c[4], c[5]), ak2^2, c[7])
		end
		@inbounds for jj = 280:min(m,2279)
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			ak2 = (BigFloat(1) / ak)^2
			Jk2[jj] = BigFloat(1) / (BigFloat(π) * ak) * muladd(ak2^2, c[5], c[7])
		end
		@inbounds for jj = 2280:m
			ak = BigFloat(π) * (jj - BigFloat(1) / BigFloat(4))
			Jk2[jj] = BigFloat(1) / (BigFloat(π) * ak) * c[7]
		end
		return Jk2
	end

    function gaulegtemp(x1::BigFloat, x2::BigFloat, n, digit)
        setprecision(digit)

        m = div(n + 1, 2)
	    xm = BigFloat(1) / BigFloat(2) * (x2 + x1)
	    xl = BigFloat(1) / BigFloat(2) * (x2 - x1)

		x = zeros(BigFloat, n)
		w = zeros(BigFloat, n)

		pp = BigFloat(0)
		@inbounds for i = 1:m
			z = cos(BigFloat(pi) * (BigFloat(i) - BigFloat(1) / BigFloat(4)) / (BigFloat(n) + BigFloat(1) / BigFloat(2)))

        	while true
				p1 = BigFloat(1)
				p2 = BigFloat(0)

				for j = 1:n
					p3 = p2
					p2 = p1
					p1 = (BigFloat(2 * j - 1) * z * p2 - BigFloat(j - 1) * p3) / BigFloat(j)
				end

				pp = BigFloat(n) * (z * p1 - p2) / (z * z - BigFloat(1))
				z1 = z
				z = z1 - p1 / pp

				if abs(z - z1) <= BigFloat(10)^(-BigFloat(3))
                	break
            	end
			end

			func = function (a)
				p1a = BigFloat(1)
				p2a = a

				for j = 1:n
					p3a = p2a
					p2a = p1a
					p1a = (BigFloat(2 * j - 1) * a * p2a - BigFloat(j - 1) * p3a) / BigFloat(j)
				end

				return p1a
			end

			z = find_zero(func, z)

			x[i] = xm - xl * z
			x[n + 1 - i] = xm + xl * z
			w[i] = BigFloat(2) * xl / ((BigFloat(1) - z * z ) * pp * pp)
			w[n + 1 - i] = w[i]
		end

        return x, w
	end

    function gauleg(func, x1::BigFloat, x2::BigFloat, n)
		x, w = gausslegendre(n)

		xm = BigFloat(1) / BigFloat(2) * (x1 + x2)
    	xr = BigFloat(1) / BigFloat(2) * (x2 - x1)

		return sum(i -> w[i] * func(xm + xr * x[i]), eachindex(x)) * xr
    end
end


text:
# Example for using Experiments

using MAB
import PyPlot

bandit  = [
    Arms.Bernoulli( 0.90 ),
    Arms.Bernoulli( 0.85 ),
    Arms.Bernoulli( 0.80 ),
    Arms.Bernoulli( 0.75 )
    # Arms.Normal( 0.36, 1.00 ),
    # Arms.Normal( 0.20, 1.00 ),
    # Arms.Normal( 0.81, 1.00 ),
    # Arms.Normal( 0.56, 1.00 ),
    # Arms.Exponential( 1.2 ),
    # Arms.Exponential( 2.3 ),
    # Arms.Exponential( 3.4 ),
    # Arms.Exponential( 4.5 )
    # Arms.Beta( 0.60, 0.40 ),
    # Arms.Sinusoidal( 800, 0 ),
    # Arms.Sinusoidal( 800, π ),
    # Arms.Sinusoidal( 800, π/2 ),
    # Arms.Sinusoidal( 800, 3*π/4 )
    # Arms.Pulse( 1000, 820, 100 ),
    # Arms.Pulse( 1000, 700, 100 ),
    # Arms.Pulse( 1000, 520, 100 )
    # Arms.Square( 200, Dict([(1,0.5),(40, 0.8),(120,0.25),(180,0.40)])),
    # Arms.Square( 200, Dict([(20, 0.3),(80,0.90),(150,0.40),(176,0.10)])),
]

noOfArms = length( bandit )

testAlgs = [
    # UniformStrategy( noOfArms ),
    # epsGreedy( noOfArms, 0.05 ),
    # epsGreedy( noOfArms, 1.00 ),
    # epsNGreedy( noOfArms, 5, 0.05 ),
    # epsNGreedy( noOfArms, 1/noOfArms, 1.0 ),
    # epsNGreedy( noOfArms ),
    # EXP3( noOfArms, 0.05 ),
    # EXP3( noOfArms, 0.10 ),
    # UCB1( noOfArms, 1/8 ),
    # UCB1( noOfArms, 1/4 ),
    # UCB1( noOfArms, 1/2 ),
    # UCB1( noOfArms, 1/√2 ),
    # UCB1( noOfArms, 1 ),
    # UCB1( noOfArms, √2 ),
    # UCB1( noOfArms, 2 ),
    # TS( noOfArms ),
    # OTS( noOfArms ),
    # DynamicTS( noOfArms, 10 ),
    # UCBNormal( noOfArms ),
    # EXP31( noOfArms ),
    # SoftMax( noOfArms, 0.008 ),
    # SoftMax( noOfArms, 0.009 ),
    # SoftMax( noOfArms, 0.010 ),
    # UniformStrategy( noOfArms ),
    # epsGreedy( noOfArms, 0.01 ),
    # epsGreedy( noOfArms, 0.05 ),
    # EXP3( noOfArms, 0.10 ),
    # REXP3( noOfArms, 0.10, 100 ),
    # UCB1( noOfArms ),
    # DUCB( noOfArms, 1.00, 0.5 ),
    # DUCB( noOfArms, 0.99 ),
    # DUCB( noOfArms, 0.99, 5.0 ),
    # EXP3( noOfArms, 0.90 ),
    # EXP3IX( noOfArms, 0.90, 0.60 )
    # TSNormal( noOfArms, 1 ),
    # TSNormal( noOfArms, 10 ),
    # TSNormal( noOfArms, 100 )
    # UCBV( noOfArms ),
    # UCBV( noOfArms, ζ = 1, c = 1),
    # UCBV( noOfArms, ζ = 2.0, c = 1),
    # BayesUCB( noOfArms ),
    # KLUCB( noOfArms )
    # TS( noOfArms ),
    # OTS( noOfArms ),
    # dTS( noOfArms ),
    # dOTS( noOfArms ),
    # dTS( noOfArms, 0.90 ),
    # dOTS( noOfArms, 0.90 ),
    # TSGaussPrior( noOfArms )
    GradientBandit( noOfArms, 0.1 ),
    KLMANB( noOfArms, 2.0, 0.0 )
]

exp1 = Experiments.Compare( bandit, testAlgs )
# run
noOfRounds      = 1000
noOfTimeSteps   = 2000
result = Experiments.run( exp1, noOfTimeSteps, noOfRounds )
fig = PyPlot.figure()
## Plot avg rewards
for alg in keys(result)
    PyPlot.plot( 1:noOfTimeSteps, cumsum(result[alg], dims=1)./collect(1:noOfTimeSteps), label = alg )
end

PyPlot.xlabel( "Timesteps" )
PyPlot.ylabel( "Avg. Reward" )
PyPlot.title( "Average Reward (Normalized for $noOfRounds rounds)")
# ax = PyPlot.gca()
# ax[:set_ylim]( [0.00,1.00] )
PyPlot.legend()
PyPlot.grid()
PyPlot.show()


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

include("delays.jl")

using Actors, Test, .Delays

t = Ref{Task}()                # this is for debugging

arg = Args(1, 2, c=3, d=4)
@test arg.args == (1, 2)
@test arg.kwargs == pairs((c=3, d=4))

a = b = c = d = 1
e = f = 2

incx(x, by; y=0, z=0) = x+y+z + by
subx(x, y, sub; z=0) = x+y+z - sub

me = Actors.newLink()
A = Actors.spawn(Bhv(incx, a, y=b, z=c), taskref=t)
sleep(0.5)
@test t[].state == :runnable

# test diag and actor startup, become! (implicitly)
act = request(A, Actors.Diag, :act)
@test act.sta == nothing
@test act.bhv.f == incx
@test act.bhv.a == (1,)
@test act.bhv.kw == pairs((y=1,z=1))

# test explicitly become!
become!(A, subx, a, b, z=c)
@test @delayed act.bhv.f == subx
@test act.bhv.a == (1,1)
@test act.bhv.kw == pairs((z=1,))

# test update!
update!(A, (1, 2, 3))
@test @delayed act.sta == (1,2,3)
update!(A, Args(2,3, x=1, y=2), s=:arg)
@test @delayed act.bhv.a == (2,3)
@test act.bhv.kw == pairs((x=1,y=2,z=1))
update!(A, :dummy, s=:mode)
@test @delayed act.mode == :dummy
@test A.mode == :dummy

# test query
query(A, me)
@test receive(me).y == (1,2,3)
@test query(A) == (1,2,3)
@test query(A, :res) == nothing
@test query(A, :bhv).f == subx

# test call
become!(A, incx, a, y=b, z=c)
@test call(A, 1) == 4
@test query(A, :res) == 4
@test query(A) == (1,2,3)

# test cast
cast(A, 2)
@test query(A, :res) == 5
update!(A, Args(5, y=1,z=1), s=:arg)
cast(A, 3)
@test query(A, :res) == 10

update!(A, Args(a, y=3,z=3))
cast(A, 3)
@test query(A, :res) == 10
@test query(A) == (1,2,3)

# test exec
exec(A, me, Bhv(cos, 2pi))
@test receive(me).y == 1
exec(A, me, sin, 2pi)
@test receive(me).y == sin(2pi)
@test exec(A, Bhv(cos, 2pi)) == 1
@test exec(A, sin, 2pi) == sin(2pi)

# test init!
init!(A, cos, 2pi)
@test @delayed act.init.f == cos

# test term!
tvar = [:ndef]
term(x) = tvar[1] = x
term!(A, term)
@test @delayed act.term.f == term

# test exit!
exit!(A)
@test @delayed t[].state == :done
@test A.chn.state == :closed
@test tvar[1] == :normal


text:
######
###### `inv`
######

function frule((_, Δx), ::typeof(inv), x::AbstractArray)
    Ω = inv(x)
    return Ω, -Ω * Δx * Ω
end

function rrule(::typeof(inv), x::AbstractArray)
    Ω = inv(x)
    function inv_pullback(ΔΩ)
        return NoTangent(), -Ω' * ΔΩ * Ω'
    end
    return Ω, inv_pullback
end

#####
##### `*`
#####

function rrule(
    ::typeof(*),
    A::AbstractVecOrMat{<:CommutativeMulNumber},
    B::AbstractVecOrMat{<:CommutativeMulNumber},
)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        return (
            NoTangent(),
            InplaceableThunk(
                @thunk(Ȳ * B'),
                X̄ -> mul!(X̄, Ȳ, B', true, true)
            ),
            InplaceableThunk(
                @thunk(A' * Ȳ),
                X̄ -> mul!(X̄, A', Ȳ, true, true)
            )
        )
    end
    return A * B, times_pullback
end

function rrule(
    ::typeof(*),
    A::AbstractVector{<:CommutativeMulNumber},
    B::AbstractMatrix{<:CommutativeMulNumber},
)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        @assert size(B, 1) === 1   # otherwise primal would have failed.
        return (
            NoTangent(),
            InplaceableThunk(
                @thunk(Ȳ * vec(B')),
                X̄ -> mul!(X̄, Ȳ, vec(B'), true, true)
            ),
            InplaceableThunk(
                @thunk(A' * Ȳ),
                X̄ -> mul!(X̄, A', Ȳ, true, true)
            )
        )
    end
    return A * B, times_pullback
end

function rrule(
   ::typeof(*), A::CommutativeMulNumber, B::AbstractArray{<:CommutativeMulNumber}
)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        return (
            NoTangent(),
            @thunk(dot(Ȳ, B)'),
            InplaceableThunk(
                @thunk(A' * Ȳ),
                X̄ -> mul!(X̄, conj(A), Ȳ, true, true)
            )
        )
    end
    return A * B, times_pullback
end

function rrule(
    ::typeof(*), B::AbstractArray{<:CommutativeMulNumber}, A::CommutativeMulNumber
)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        return (
            NoTangent(),
            InplaceableThunk(
                @thunk(A' * Ȳ),
                X̄ -> mul!(X̄, conj(A), Ȳ, true, true)
            ),
            @thunk(dot(Ȳ, B)'),
        )
    end
    return A * B, times_pullback
end

#####
##### `muladd`
#####

function rrule(
        ::typeof(muladd),
        A::AbstractMatrix{<:CommutativeMulNumber},
        B::AbstractVecOrMat{<:CommutativeMulNumber},
        z::Union{CommutativeMulNumber, AbstractVecOrMat{<:CommutativeMulNumber}},
    )
    # The useful case, mul! fused with +
    function muladd_pullback_1(ȳ)
        Ȳ = unthunk(ȳ)
        matmul = (
            InplaceableThunk(
                @thunk(Ȳ * B'),
                dA -> mul!(dA, Ȳ, B', true, true)
            ),
            InplaceableThunk(
                @thunk(A' * Ȳ),
                dB -> mul!(dB, A', Ȳ, true, true)
            )
        )
        addon = if z isa Bool
            NoTangent()
        elseif z isa Number
            @thunk(sum(Ȳ))
        else
            InplaceableThunk(
                @thunk(sum!(similar(z, eltype(Ȳ)), Ȳ)),
                dz -> sum!(dz, Ȳ; init=false)
            )
        end
        (NoTangent(), matmul..., addon)
    end
    return muladd(A, B, z), muladd_pullback_1
end

function rrule(
        ::typeof(muladd),
        ut::LinearAlgebra.AdjOrTransAbsVec{<:CommutativeMulNumber},
        v::AbstractVector{<:CommutativeMulNumber},
        z::CommutativeMulNumber,
    )
    # This case is dot(u,v)+z, but would also match signature above.
    function muladd_pullback_2(ȳ)
        dy = unthunk(ȳ)
        ut_thunk = InplaceableThunk(
            @thunk(v' .* dy),
            dut -> dut .+= v' .* dy
        )
        v_thunk = InplaceableThunk(
            @thunk(ut' .* dy),
            dv -> dv .+= ut' .* dy
        )
        (NoTangent(), ut_thunk, v_thunk, z isa Bool ? NoTangent() : dy)
    end
    return muladd(ut, v, z), muladd_pullback_2
end

function rrule(
        ::typeof(muladd),
        u::AbstractVector{<:CommutativeMulNumber},
        vt::LinearAlgebra.AdjOrTransAbsVec{<:CommutativeMulNumber},
        z::Union{CommutativeMulNumber, AbstractVecOrMat{<:CommutativeMulNumber}},
    )
    # Outer product, just broadcasting
    function muladd_pullback_3(ȳ)
        Ȳ = unthunk(ȳ)
        proj = (
            @thunk(vec(sum(Ȳ .* conj.(vt), dims=2))),
            @thunk(vec(sum(u .* conj.(Ȳ), dims=1))'),
        )
        addon = if z isa Bool
            NoTangent()
        elseif z isa Number
            @thunk(sum(Ȳ))
        else
            InplaceableThunk(
                @thunk(sum!(similar(z, eltype(Ȳ)), Ȳ)),
                dz -> sum!(dz, Ȳ; init=false)
            )
        end
        (NoTangent(), proj..., addon)
    end
    return muladd(u, vt, z), muladd_pullback_3
end

#####
##### `/`
#####

function rrule(::typeof(/), A::AbstractVecOrMat{<:Real}, B::AbstractVecOrMat{<:Real})
    Aᵀ, dA_pb = rrule(adjoint, A)
    Bᵀ, dB_pb = rrule(adjoint, B)
    Cᵀ, dS_pb = rrule(\, Bᵀ, Aᵀ)
    C, dC_pb = rrule(adjoint, Cᵀ)
    function slash_pullback(Ȳ)
        # Optimization note: dAᵀ, dBᵀ, dC are calculated no matter which partial you want
        _, dC = dC_pb(Ȳ)
        _, dBᵀ, dAᵀ = dS_pb(unthunk(dC))

        ∂A = last(dA_pb(unthunk(dAᵀ)))
        ∂B = last(dA_pb(unthunk(dBᵀ)))

        (NoTangent(), ∂A, ∂B)
    end
    return C, slash_pullback
end

#####
##### `\`
#####

function rrule(::typeof(\), A::AbstractVecOrMat{<:Real}, B::AbstractVecOrMat{<:Real})
    Y = A \ B
    function backslash_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        ∂A = @thunk begin
            B̄ = A' \ Ȳ
            Ā = -B̄ * Y'
            Ā = add!!(Ā, (B - A * Y) * B̄' / A')
            Ā = add!!(Ā, A' \ Y * (Ȳ' - B̄'A))
            Ā
        end
        ∂B = @thunk A' \ Ȳ
        return NoTangent(), ∂A, ∂B
    end
    return Y, backslash_pullback

end

#####
##### `\`, `/` matrix-scalar_rule
#####

function rrule(::typeof(/), A::AbstractArray{<:CommutativeMulNumber}, b::CommutativeMulNumber)
    Y = A/b
    function slash_pullback_scalar(ȳ)
        Ȳ = unthunk(ȳ)
        Athunk = InplaceableThunk(
            @thunk(Ȳ / conj(b)),
            dA -> dA .+= Ȳ ./ conj(b),
        )
        bthunk = @thunk(-dot(A,Ȳ) / conj(b^2))
        return (NoTangent(), Athunk, bthunk)
    end
    return Y, slash_pullback_scalar
end

function rrule(::typeof(\), b::CommutativeMulNumber, A::AbstractArray{<:CommutativeMulNumber})
    Y, back = rrule(/, A, b)
    function backslash_pullback(dY)  # just reverses the arguments!
        d0, dA, db = back(dY)
        return (d0, db, dA)
    end
    return Y, backslash_pullback
end

#####
##### Negation (Unary -)
#####

function rrule(::typeof(-), x::AbstractArray)
    function negation_pullback(ȳ)
        return NoTangent(), InplaceableThunk(@thunk(-ȳ), ā -> ā .-= ȳ)
    end
    return -x, negation_pullback
end

#####
##### Addition (Multiarg `+`)
#####

function rrule(::typeof(+), arrs::AbstractArray...)
    y = +(arrs...)
    arr_axs = map(axes, arrs)
    function add_pullback(dy)
        return (NoTangent(), map(ax -> reshape(dy, ax), arr_axs)...)
    end
    return y, add_pullback
end


text:
using BenchmarkTools

flat(arr) = mapreduce(x -> x == [] || x[1] === x ? x : flat(x), vcat, arr, init=[])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
##
reload("MetabolicFluxAnalysis")
using MetabolicFluxAnalysis

# load file
fluxi_1_f = open("datasets/antoniewicz_toy_network/toy_network_01.fluxi")
fluxi_1_s = readstring(fluxi_1_f)
fluxi1 = MetabolicFluxAnalysis.parseFluxi( fluxi_1_s )

# create an evaluator for the forward simulation
ilee_1 = MetabolicFluxAnalysis.create_julia_ile_evaluator( fluxi1.net , fluxi1.substrate_configs["SC_010"] ,fluxi1.measurements["MassDistribution_01"] )

# create flux polytope:

fp1 = MetabolicFluxAnalysis.get_flux_polyhedron( fluxi1.net , fluxi1.linear_constraints["ConstraintsExchange"] , max_val=200 )

# parametrize the flux polytope
rfd_parametrization = MetabolicFluxAnalysis.get_rounded_fulldimensional_parametrization( fp1 ; remove_red_constr=1)

@assert MetabolicFluxAnalysis.isinside( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt[:,:] )

#                    v1  v2     v3   v4   v5   v6      Ain  Eout  Fout
reference_flux = [  100.; 110 ; 50 ; 20 ; 20 ; 80   ;  100 ; 60 ; 80 ]

ilee_1.f_fwsim( ilee_1.substrate_mds , reference_flux[:,[1;1;1;1]])

md_xx  = MetabolicFluxAnalysis.eval_flux2mds( ilee_1 , repmat(reference_flux[:,:],1,100) )
sse_xx = MetabolicFluxAnalysis.eval_flux2sse( ilee_1 , repmat(reference_flux[:,:],1,100) )

# generate some unifrom samples over polytope:
xx = MetabolicFluxAnalysis.sample_unif_hr( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt , 5000 , 1000 );
xx_flux = MetabolicFluxAnalysis.eval_opt2raw(rfd_parametrization,xx)

sse_xx = MetabolicFluxAnalysis.eval_flux2sse( ilee_1 , xx_flux )
#sse_xx = MetabolicFluxAnalysis.eval_flux2sse( ilee_1 , reference_flux[:,:] )

using Plots
gr()
Plots.scatter( xx_flux[2,:],xx_flux[5,:] , zcolor=sse_xx[:],c=:lightrainbow,ms=4.0,marker=(stroke(0,:gray)),clims=(0,20))

#Plots.plot!()

using Plots
gr()
plot(xx_flux[2,:],xx_flux[5,:],sse_xx,nlevels=10)
ipX = ipY = linspace(-3, 3, 7)
ipZ = GR.interp2(xx_flux[2,:],xx_flux[5,:],sse_xx,ipX,ipY )

default(size=(500,300), leg=true)
y = rand(100)
scatter(y, z=abs(y-.5), m=10, c=:heat, lab="grad")
plot!(0:10:100,rand(11,4),lab="lines", palette=:grays)

# benchmark of unif. sampling:
tic(); xx = MetabolicFluxAnalysis.sample_unif_hr( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt , 1000 ,1000 ); toc();

using MetabolicFluxAnalysisWorkflows
f_cost = (x) -> MetabolicFluxAnalysis.eval_flux2sse( ilee_1 , MetabolicFluxAnalysis.eval_opt2raw(rfd_parametrization,x[:,:]) )[1]
opti_result = MetabolicFluxAnalysisWorkflows.find_optimum( rfd_parametrization.p_opt  , f_cost  , xx[:,1] )
opti_flux = MetabolicFluxAnalysis.eval_opt2raw( rfd_parametrization , opti_result[2] )
scatter!(opti_flux[2:2],opti_flux[5:5])
opti_result = MetabolicFluxAnalysisWorkflows.find_optimum( rfd_parametrization.p_opt  , f_cost  , xx )

#julia_fwsim_1 =  MetabolicFluxAnalysis.compile_julia_fwsim( fluxi1.net , fluxi1.measurements["MassDistribution_01"] )
#emudecomp_1 = MetabolicFluxAnalysis.compute_emu_decomp( fluxi1.net , [ fluxi1.measurements["MassDistribution_01"][1].sass ] )
#emu_graph_1 = MetabolicFluxAnalysis.create_emu_graph(emudecomp_1)
#fluxi_ast_1 = MetabolicFluxAnalysis.create_ast_from_emu_graph( emu_graph_1 )
#julia_fw_sim_1 = MetabolicFluxAnalysis.ast_to_juliacode_01(fluxi1,fluxi_ast_1,target_smrs)

#str_json_01 = MetabolicFluxAnalysis.export_json(fluxi1)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "pciutils"
version = v"3.7.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/pciutils/pciutils.git", "864aecdea9c7db626856d8d452f6c784316a878c")
]

dependencies = Dependency[
]

# Bash recipe for building across all platforms
script = raw"""
cd pciutils
make install PREFIX=${prefix} SHARED=yes SBINDIR=${bindir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> !Sys.iswindows(p) && !Sys.isapple(p), supported_platforms(;experimental=true))

# The products that we will ensure are always built
products = Product[
    ExecutableProduct("lspci", :lspci),
    ExecutableProduct("setpci", :setpci),
    LibraryProduct("libpci", :libpci)
]

# Build the tarballs, and possibly a `build.jl` as well.
# gcc7 constraint from boost
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using DataFrames
using JointMoments

include("data.jl")

const ITERMAX = 25

function timing()
    time_cov = true
    time_coskew = true
    time_cokurt = true
    datasets = (data, bigdata)
    # datasets = (bigdata,)
    for (idx, dataset) in enumerate(datasets)
        df = DataFrame()
        rows, cols = size(dataset)
        label = string("Data: ", rows, "x", cols, " ", typeof(dataset), "\n")
        print_with_color(:white, label)
        for i = 1:2
            itermax = (i == 1) ? 1 : ITERMAX
            row = {
                :bias => 1,
                :flatten => false,
                :standardize => true,
                :function => :_cov,
                :elapsed => 0.0,
                :elapsed_error => 0.0,
            }
            # parameters: dense, standardize, flatten
            params = [
                (true, true, true),
                (true, true, false),
                (true, false, true),
                (false, true, true),
                (true, false, false),
                (false, true, false),
                (false, false, true),
                (false, false, false),
            ]
            for p in params
                row[:dense] = p[1]
                if time_cov
                    row[:standardize] = NaN
                    row[:flatten] = false
                    row[:function] = :_cov
                    elapse = (Float64)[]
                    for j = 1:itermax
                        push!(elapse,
                            @elapsed _cov(
                                dataset,
                                dense=row[:dense],
                                bias=row[:bias],
                            )
                        )
                    end
                    row[:elapsed] = mean(elapse)
                    row[:elapsed_error] = std(elapse) / ITERMAX
                    if i == 2
                        if size(df) == (0,0)
                            for (key, value) in row
                                df[key] = value
                            end
                        else
                            push!(df, [value for (key, value) in row])
                        end
                    end
                end
                if time_coskew
                    row[:dense] = p[1]
                    row[:standardize] = p[2]
                    row[:flatten] = p[3]
                    row[:function] = :coskew
                    elapse = (Float64)[]
                    for j = 1:itermax
                        push!(elapse,
                            @elapsed coskew(
                                dataset,
                                standardize=row[:standardize],
                                flatten=row[:flatten],
                                bias=row[:bias],
                                dense=row[:dense],
                            )
                        )
                    end
                    row[:elapsed] = mean(elapse)
                    row[:elapsed_error] = std(elapse) / ITERMAX
                    if i == 2
                        if size(df) == (0,0)
                            for (key, value) in row
                                df[key] = value
                            end
                        else
                            push!(df, [value for (key, value) in row])
                        end
                    end
                end
                if time_cokurt
                    row[:dense] = p[1]
                    row[:standardize] = p[2]
                    row[:flatten] = p[3]
                    row[:function] = :cokurt
                    elapse = (Float64)[]
                    for j = 1:itermax
                        push!(elapse,
                            @elapsed cokurt(
                                dataset,
                                standardize=row[:standardize],
                                flatten=row[:flatten],
                                bias=row[:bias],
                                dense=row[:dense],
                            )
                        )
                    end
                    row[:elapsed] = mean(elapse)
                    row[:elapsed_error] = std(elapse) / ITERMAX
                    if i == 2
                        if size(df) == (0,0)
                            for (key, value) in row
                                df[key] = value
                            end
                        else
                            push!(df, [value for (key, value) in row])
                        end
                    end
                end
            end
        end
        df = sort(df, cols=:elapsed)
        println(df)
    end
end

timing()


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ a9c917f0-37ec-11eb-39c1-8d9fec2fdaf1
# can take a while....
# might need to install DifferentialEquations & Plots with
# ]add Differential Equations Plots
begin
	using DifferentialEquations
	using Plots
	using Base.Iterators: product
end

# ╔═╡ 7866daf0-5eaa-11eb-2c25-afec8b81fb56
md"
## Investment problem

```math
\max_{i(t)} \int_0^T e^{-\rho t} \left\{p k(t) - c i(t) - \frac{d}{2} i(t)^2 \right\} dt
\qquad st \qquad
\dot k(t) = i(t) - \delta k(t)
```
given boundary conditions
```math
k(0) = k_0 \qquad k(T) = k_T \qquad T:\text{ given}
```

Necessary conditions
```math
\begin{align*}
\dot k &= i - \delta k \\
\dot \psi &= (\rho + \delta)\psi - p
\end{align*}
```
"

# ╔═╡ 6c39a320-5eaa-11eb-0729-0fdeb7f66e7f
# tell plots to use GR b/c it's fastest
gr();

# ╔═╡ bcf87142-5ea7-11eb-2359-115b8572fbe2
begin
	# parameters
	px=0.5
	c=3
	d=1
	ρ=0.07
	δ=0.05

	# steady state
	ssi = (px / (ρ+δ) - c) / d
	ssk = ssi/δ
end

# ╔═╡ 6b845912-37eb-11eb-3948-d513425270f2
begin
	doti(k,i) = (ρ+δ)*(i - ssi)
	dotk(k,i) = i - δ*k

	# solvers in DifferentialEquations.jl require functions that
	# take the form   `f!(dy, y, params, t)`
	function dot_ki!(dy, y, parm, t)
		k,i = y
		dy[1] = dotk(k,i)
		dy[2] = doti(k,i)
		return dy
	end
end

# ╔═╡ 25be1680-37f2-11eb-2c06-99283062577c
# create a "mesh"
begin
	# grid of points in K, I space
	kspace = range(ssk/2, stop = 1.5*ssk, length=15)
	ispace = range(ssi/2, stop = 1.5*ssi, length=15)

	kispace = product(kspace, ispace)
	KK = [k for (k,i) in kispace]
	II = [i for (k,i) in kispace]

	SCALE = 2
	dotII = [doti(k,i) for (k,i) in kispace].*SCALE
	dotKK = [dotk(k,i) for (k,i) in kispace].*SCALE
end

# ╔═╡ 001a0000-5ea9-11eb-1a57-29544ffeb338
begin
	plt_quiv = plot(;legend=:bottomright, xlabel="\$k\$", ylabel="\$i\$")

	# motion vectors
	# note, use `vec` b/c need to input a vector
	quiver!(plt_quiv, vec(KK), vec(II); quiver=(vec(dotKK),vec(dotII)), alpha=0.2)

	# steady state
	# also note that we have to give vectors to scatter
	scatter!(plt_quiv, [ssk],[ssi], label="SS")

	# nullclines
	plot!(plt_quiv, kspace, k -> k*δ, label="\$\\dot k = 0\$")
	plot!(plt_quiv, kspace, k -> ssi, label="\$\\dot i = 0\$")
end

# ╔═╡ d412fe80-37f2-11eb-3c13-87530a5aa318
begin
	# initial conditions
	# (k₀ i₀)
	y0 = [12.0, 1.19]

	# time-span for problem
	# NOTE: this MUST be in Floats, not Ints!
	tmax = 25.0
	tspan = (0.0, tmax)

	# define a problem
	prob = ODEProblem(dot_ki!, y0, tspan);
end

# ╔═╡ e9758ae2-37f2-11eb-2299-c17952d5c072
# solve the problem using the `Tsit5()` method
# Note how much better our tolerances are
sol = solve(prob, Tsit5())

# ╔═╡ 587a34e0-5ea9-11eb-09ee-b9c836837819
begin
	TT = 0 : 0.1 : tmax
	kpath = first.(sol.(TT))
	ipath = last.(sol.(TT))

	# need to copy plot b/c of updating in Pluto
	pltq = deepcopy(plt_quiv)
	plot!(pltq, kpath, ipath; label="\$(k_t, i_t)\$")
end

# ╔═╡ faaf5290-5ea9-11eb-3761-7db74d6930e8
plot(
    plot(TT, t -> sol(t)[1], title="\$k(t)\$", legend=:none),
    plot(TT, t -> sol(t)[2], title="\$i(t)\$", legend=:none),
    plot(
        sol,
        labels=["\$k(t)\$" "\$i(t)\$"], # NOTE: labels are row matrix!
        legend=:right,
        title="plot(sol)"
    ),
    layout=(1,3)
)

# ╔═╡ Cell order:
# ╟─7866daf0-5eaa-11eb-2c25-afec8b81fb56
# ╠═a9c917f0-37ec-11eb-39c1-8d9fec2fdaf1
# ╠═6c39a320-5eaa-11eb-0729-0fdeb7f66e7f
# ╠═bcf87142-5ea7-11eb-2359-115b8572fbe2
# ╠═6b845912-37eb-11eb-3948-d513425270f2
# ╠═25be1680-37f2-11eb-2c06-99283062577c
# ╠═001a0000-5ea9-11eb-1a57-29544ffeb338
# ╠═d412fe80-37f2-11eb-3c13-87530a5aa318
# ╠═e9758ae2-37f2-11eb-2299-c17952d5c072
# ╠═587a34e0-5ea9-11eb-09ee-b9c836837819
# ╠═faaf5290-5ea9-11eb-3761-7db74d6930e8


text:
import Base: |, merge
export ←, |

################################## Exact Conditioning ######################################

"""
    Observation

Represents fixing a paricular (finite) GP to have a particular (vector) value.
"""
struct Observation{Tf<:AbstractGP, Ty<:AbstractVector}
    f::Tf
    y::Ty
end
←(f, y) = Observation(f, y)
get_f(c::Observation) = c.f
get_y(c::Observation) = c.y
function merge(c::Union{AbstractVector{<:Observation}, Tuple{Vararg{Observation}}})
    return BlockGP([get_f.(c)...])←BlockVector([get_y.(c)...])
end

"""
    |(g::AbstractGP, c::Observation)

Condition `g` on observation `c`.
"""
|(g::GP, c::Observation) = GP(|, g, c.f, CondCache(mean_vec(c.f), cov(c.f), c.y))
|(g::BlockGP, c::Observation) = BlockGP(g.fs .| Ref(c))

function μ_p′(::typeof(|), g::AbstractGP, f::AbstractGP, cache::CondCache)
    return iszero(kernel(f, g)) ?
        mean(g) :
        ConditionalMean(cache, mean(g), kernel(f, g))
end
function k_p′(::typeof(|), g::AbstractGP, f::AbstractGP, cache::CondCache)
    return iszero(kernel(f, g)) ?
        kernel(g) :
        ConditionalKernel(cache, kernel(f, g), kernel(g))
end
function k_p′p(::typeof(|), g::AbstractGP, f::AbstractGP, cache::CondCache, h::AbstractGP)
    return iszero(kernel(f, g)) || iszero(kernel(f, h)) ?
        kernel(g, h) :
        ConditionalCrossKernel(cache, kernel(f, g), kernel(f, h), kernel(g, h))
end
function k_pp′(h::AbstractGP, ::typeof(|), g::AbstractGP, f::AbstractGP, cache::CondCache)
    return iszero(kernel(f, g)) || iszero(kernel(f, h)) ?
        kernel(h, g) :
        ConditionalCrossKernel(cache, kernel(f, h), kernel(f, g), kernel(h, g))
end

# Sugar
|(g::AbstractGP, c::Tuple{Vararg{Observation}}) = g | merge(c)
|(g::Tuple{Vararg{AbstractGP}}, c::Observation) = deconstruct(BlockGP([g...]) | c)
function |(g::Tuple{Vararg{AbstractGP}}, c::Tuple{Vararg{Observation}})
    return deconstruct(BlockGP([g...]) | merge(c))
end

################################# Approximate Conditioning #################################

abstract type AbstractConditioner end

FillArrays.Zeros(x::BlockVector) = BlockVector(Zeros.(x.blocks))

"""
    Titsias <: AbstractConditioner

Construct an object which is able to compute an approximate posterior.
"""
struct Titsias{Tu<:AbstractGP, Tm<:AV{<:Real}, Tγ} <: AbstractConditioner
    u::Tu
    m′u::Tm
    γ::Tγ
    function Titsias(u::Tu, m′u::Tm, Σ′uu::AM) where {Tu, Tm}
        μ = EmpiricalMean(Zeros(m′u))
        Σ = EmpiricalKernel(Xtinv_A_Xinv(Σ′uu, cov(u)))
        γ = GP(μ, Σ, u.gpc)
        return new{Tu, Tm, typeof(γ)}(u, m′u, γ)
    end
end
function |(g::GP, c::Titsias)
    g′ = g | (c.u←c.m′u)
    ĝ = project(kernel(c.u, g), c.γ)
    return g′ + ĝ
end

function optimal_q(f::AbstractGP, y::AV{<:Real}, u::AbstractGP, σ::Real)
    μᵤ, Σᵤᵤ = mean_vec(u), cov(u)
    U = cholesky(Σᵤᵤ).U
    Γ = broadcast(/, U' \ xcov(u, f), σ)
    Ω, δ = LazyPDMat(Symmetric(Γ * Γ' + I), 0), y - mean_vec(f)
    Σ′ᵤᵤ = Xt_invA_X(Ω, U)
    μ′ᵤ = μᵤ + broadcast(/, U' * (Ω \ (Γ * δ)), σ)
    return μ′ᵤ, Σ′ᵤᵤ
end
optimal_q(c::Observation, u::AbstractGP, σ::Real) = optimal_q(c.f, c.y, u, σ)

# Sugar.
function optimal_q(f::AV{<:AbstractGP}, y::AV{<:AV{<:Real}}, u::AbstractGP, σ::Real)
    return optimal_q(BlockGP(f), BlockVector(y), u, σ)
end
function optimal_q(f::AbstractGP, y::AV{<:Real}, u::AV{<:AbstractGP}, σ::Real)
    return optimal_q(f, y, BlockGP(u), σ)
end
function optimal_q(f::AV{<:AbstractGP}, y::AV{<:AV{<:Real}}, u::AV{<:AbstractGP}, σ::Real)
    return optimal_q(BlockGP(f), BlockVector(y), BlockGP(u), σ)
end

|(g::Tuple{Vararg{AbstractGP}}, c::Titsias) = deconstruct(BlockGP([g...]) | c)
function |(g::BlockGP, c::Titsias)
    return BlockGP(g.fs .| Ref(c))
end

"""
    Titsias(
        c::Union{Observation, Vector{<:Observation}},
        u::Union{AbstractGP, Vector{<:AbstractGP}},
        σ::Real,
    )

Instantiate the saturated Titsias conditioner.
"""
function Titsias(c::Observation, u::AbstractGP, σ::Real)
    return Titsias(u, optimal_q(c, u, σ)...)
end
function Titsias(c::Vector{<:Observation}, u::Vector{<:AbstractGP}, σ::Real)
    return Titsias(merge(c), BlockGP(u), σ)
end
function Titsias(c::Observation, u::Vector{<:AbstractGP}, σ::Real)
    return Titsias(c, BlockGP(u), σ::Real)
end
function Titsias(c::Vector{<:Observation}, u::AbstractGP, σ::Real)
    return Titsias(merge(c), u,  σ)
end


text:
# Analyze convergence of Monte Carlo tree search for various horizon lenghts

using SubmodularMaximization
using PyPlot
using Statistics
using Base.Threads
using Base.Iterators
using JLD2
using Printf

close("all")

println("Loading rollout library")
library_file = "./data/rollout_library.jld2"
@time @load library_file data num_robots trials
rollouts = data

# Create a map from number of robots to a vector of configs sample pairs
samples = Dict(map(num_robots) do num_robots
                 num_robots => vcat(map(trials) do trial
                                      rollout = rollouts[num_robots, trial]

                                      [(configs=rollout.configs, sample=sample)
                                       for sample in rollout.samples]
                                    end...)
               end)

# Everything has the same number of samples
sample_inds = 1:length(first(values(samples)))

data_file = "./data/multi_robot_mcts_convergence_data.jld2"
reprocess = false

num_mcts_samples = 2 .^ (4:13)
horizons = 1:6

solver = solve_sequential

all_tests = product(horizons, num_mcts_samples, num_robots, sample_inds)
all_configurations = product(horizons, num_mcts_samples, num_robots)

function get_data()

  # Both guarded by the mutex
  results_mutex = Mutex()
  results = Dict()
  num_tests_completed = 0

  if !isfile(data_file) || reprocess

    iters = collect(all_tests)
    @threads for (horizon, num_mcts_samples, num_robots,
                  sample_ind) in iters

      test = samples[num_robots][sample_ind]

      sample = test.sample

      configs = MultiRobotTargetTrackingConfigs(test.configs;
                                          horizon = horizon,
                                          solver_iterations = num_mcts_samples)

      problem = MultiRobotTargetTrackingProblem(sample.robot_states,
                                                sample.target_filters,
                                                configs)
      solution = solver(problem)
      reward = solution.value

      lock(results_mutex)

      results[horizon, num_mcts_samples, num_robots, sample_ind] = reward
      num_tests_completed +=1

      unlock(results_mutex)

      println("Horizon: ", horizon,
              " Num. samples: ", num_mcts_samples,
              " Num. robots: ", num_robots,
              " Trial: ", sample_ind,
              " (done, ", num_tests_completed, "/", length(all_tests),
              @sprintf(" %0.2f", 100num_tests_completed/length(all_tests)), "%)"
             )
    end

    @save data_file results
  else
    @load data_file results
  end

  results
end

println("Processing rollouts")
results = get_data()

println("Analyzing results")

#
# normalize rewards
#

# Compute the max reward over horizons/num_robots/sample_inds
max_rewards = Dict(map(product(horizons, num_robots, sample_inds)) do x
                     horizon, num_robots, sample_ind = x

                     x => maximum(num_mcts_samples) do mcts_samples
                       results[horizon, mcts_samples, num_robots, sample_ind]
                     end
                   end)

# Normalize rewards over results for different numbers of samples
normalized_rewards = Dict(map(all_tests) do x
                            horizon,num_mcts_samples,num_robots,sample_ind = x
                            configuration = (horizon, num_mcts_samples,
                                             num_robots)

                            max_reward = max_rewards[horizon, num_robots,
                                                     sample_ind]

                            x=>results[x...]/max_reward
                          end)

# Mean reward over all samples for each configuration
mean_rewards = Dict(map(all_configurations) do x
                      x => mean(ind->normalized_rewards[x..., ind], sample_inds)
                    end)

for num_robots in num_robots
  figure()
  for horizon in horizons
    horizon_data = map(num_mcts_samples) do mcts_samples
      mean_rewards[horizon, mcts_samples, num_robots]
    end

    # Print the maximum reward for that horizon and such
    println("Horizon: ", horizon, " Max reward: ", maximum(horizon_data))

    label = string("Horz.: ", horizon)
    loglog(num_mcts_samples, horizon_data, label=label)
  end

  title(string("MCTS convergence (", num_robots, " robots)"))
  xlabel("Num rollouts")
  ylabel("Frac. of max. reward")
  legend(loc="upper left")

  fig_name =  string("mcts_convergence_", num_robots, "_robots")
  save_fig("fig", fig_name)
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
#! /usr/bin/julia

# Rosetta Code, Sorting algorithms/Pancake sort

function pancakesort!{T<:Real}(a::Array{T,1})
    len = length(a)
    if len < 2
        return a
    end
    for i in len:-1:2
        j = indmax(a[1:i])
        i != j || continue
        a[1:j] = reverse(a[1:j])
        a[1:i] = reverse(a[1:i])
    end
    return a
end

pancakesort{T<:Real}(a::Array{T,1}) = pancakesort!(copy(a))

println("Testing a pancake sort.")
a = [rand(-100:100) for i in 1:20]
println("Pre:   ", a)
b = pancakesort(a)
println("Post:  ", b)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
import Elemental
const El = Elemental

using MPI

n0 = 50
n1 = 50
testNative = true
display = false
worldRank = MPI.Comm_rank(MPI.COMM_WORLD)
worldSize = MPI.Comm_size(MPI.COMM_WORLD)
if worldRank == 0
  print("worldSize=$worldSize\n")
end

using DistributedArrays

function randDArray(m,n,sparsity=0.01)
    DA = DistributedArrays.DArray((m,n)) do I
        sprandn(length(I[1]), length(I[2]), sparsity)
    end
    A = El.DistSparseMatrix{Float64}(DA)
    return A
end

function stackedFD2D(n0, n1)
    height = 2*n0*n1
    width = n0*n1
    A = El.DistSparseMatrix(Float64, height, width)
    localHeight = El.localHeight(A)
    El.reserve(A, 6*localHeight)

    for sLoc in 1:localHeight
        s = El.globalRow(A, sLoc)
        if s <= n0*n1
            x0 = ((s-1) % n0) + 1
            x1 = div((s-1), n0) + 1
            El.queueUpdate(A, s, s, 11.0)
            if x0 > 1
                El.queueUpdate(A, s, s - 1, -10.0)
            end
            if x0 < n0
                El.queueUpdate(A, s, s + 1, 20.0)
            end
            if x1 > 1
                El.queueUpdate(A, s, s - n0, -30.0)
            end
            if x1 < n1
                El.queueUpdate(A, s, s + n0, 40.0)
            end
        else
            sRel = s - n0*n1
            x0 = ((sRel-1) % n0) + 1
            x1 = div(sRel-1, n0) + 1
            El.queueUpdate(A, s, sRel, -20.0)
            if x0 > 1
                El.queueUpdate(A, s, sRel - 1, -1.0)
            end
            if x0 < n0
                El.queueUpdate(A, s, sRel + 1, -2.0)
            end
            if x1 > 1
                El.queueUpdate(A, s, sRel - n0, -3.0)
            end
            if x1 < n1
                El.queueUpdate(A, s, sRel + n0, 3.0)
            end
        end

        # The dense last column
        El.queueUpdate(A, s, width, -div(10.0, height))
    end
    El.processQueues(A)
    return A
end

if testNative
    A = stackedFD2D(n0, n1)
else
    A = randDArray(2*n0*n1,n0*n1)
end

b = El.DistMultiVec(Float64)
El.gaussian!(b, 2*n0*n1, 1)

AHeight = El.height(A)
AWidth = El.width(A)
bHeight = El.height(b)
bWidth = El.width(b)

if display
    show(IO, A)
end
ctrl = El.LPAffineCtrl(Float64,
            mehrotraCtrl = El.MehrotraCtrl(Float64,
                            solveCtrl = El.RegSolveCtrl(Float64, progress = true),
                            print = true,
                            outerEquil = true,
                            time = true))

elapsedLAV = @elapsed x = El.lav(A, b)
elapsedLAV = @elapsed x = El.lav(A, b, ctrl)

if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    println("LAV time: $elapsedLAV seconds")
end
bTwoNorm = El.nrm2(b)
bInfNorm = El.maxNorm(b)

r = copy(b)
A_mul_B!(-1.0, A, x, 1.0, r)

rTwoNorm = El.nrm2(r)
rOneNorm = El.entrywiseNorm(r, 1)

if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    println("|| b ||_2       = $bTwoNorm")
    println("|| b ||_oo      = $bInfNorm")
    println("|| A x - b ||_2 = $rTwoNorm")
    println("|| A x - b ||_1 = $rOneNorm")
end

elapsedLS = @elapsed xLS = El.leastSquares(A, b)

if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    println("LS time: $elapsedLAV seconds")
end

rLS = copy(b)
A_mul_B!(-1.0, A, xLS, 1., rLS)
if display
    El.Display( rLS, "A x_{LS} - b" )
end

rLSTwoNorm = El.nrm2(rLS)
rLSOneNorm = El.entrywiseNorm(rLS, 1)
if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    println("|| A x_{LS} - b ||_2 = $rLSTwoNorm")
    println("|| A x_{LS} - b ||_1 = $rLSOneNorm")
end

# Require the user to press a button before the figures are closed
# commSize = El.mpi.Size( El.mpi.COMM_WORLD() )
El.Finalize()


text:
"""
Presolve is a key to accelerate the linear programming solution.

TODO: Need to check and implement the algorithms in
Andersen and Andersen. Presolving in linear programming. 1995
"""
function presolve(lp::MatOI.LPForm{T, AT, VT}) where {T, AT, VT}
    @timeit TO "row reduction" begin
        tmp_A = lp.A
        to_keep = collect(1:nrows(lp))
        # Remove empty rows
        non_empty = get_empty_rows(tmp_A)
        if length(non_empty) < nrows(lp)
            println("Presolve removed $(nrows(lp)-length(non_empty)) empty rows.")
            tmp_A = tmp_A[non_empty, :]
            to_keep = non_empty
        end
        # Remove dependent rows
        dep_rows = get_dependent_rows(tmp_A)
        if length(dep_rows) > 0
            indeprows = [i for i=1:size(tmp_A, 1) if !in(i,dep_rows)]
            println("Presolve reduced $(length(dep_rows)) rows.")
            to_keep = to_keep[indeprows]
            tmp_A = tmp_A[indeprows, :]
        end

        # @show size(lp.A)
    end
    return MatOI.LPForm{T, AT, VT}(
        lp.direction,
        lp.c,
        sparse(tmp_A),
        lp.c_lb[to_keep],
        lp.c_ub[to_keep],
        lp.v_lb,
        lp.v_ub,
    )
end

function get_empty_rows(A::AbstractArray{T, 2}) where T
    nonempty = Int[]
    m = size(A, 1)
    sizehint!(nonempty, m)

    for i = 1:m
        if norm(A[i,:]) > 1e-10
            push!(nonempty, i)
        end
    end

    return nonempty
end

"""
Eliminate linearly dependent rows

The implementation is based on
Chvatal. Linear programming. 1981

The computational efficiency can be improved by
Andersen. Finding all linearly dependent rows in large-scale linear programming. 1995
"""
function get_dependent_rows(A::AbstractArray{T, 2}) where T
    m, n = size(A)

    # Append identity matrix to the last
    A2 = [A I]

    basis = collect((1+n):(m+n))
    # e_i = sparsevec([1],[1.0],m)
    e_i = Array{Float64}(undef, m)
    pi = Array{Float64}(undef, m)

    depind = Int[]
    sizehint!(depind, m)

    for i in 1:m
        # initialize compute components
        # e_i.nzind[1] = i
        fill!(e_i, 0)
        e_i[i] = 1

        # compute pi
        pi .= A2'[basis,:] \ e_i

        is_dependent = true
        for j = 1:n
            # @show A2[:,j]
            if !in(j, basis)
                view_A2_j = @view A2[:,j]
                if abs(view_A2_j' * pi) > 1e-7
                    basis[i] = j
                    # println("Pivot: row $i column $j")
                    is_dependent = false
                    break
                end
            end
        end
        if is_dependent
            # println("Row $i is dependent.")
            push!(depind, i)
        end
    end
    # @show basis[basis .> n]
    # @show basis[basis .> n] .- n

    return depind
end

# function remove_dependent_rows!(lp::StandardLpData{CuArray})
#     m, n = size(lp.A)

#     # Append identity matrix to the last
#     A2 = [lp.A CuMatrix(Matrix(I,m,m))]

#     basis = collect((1+n):(m+n))
#     e_i = CuArray{Float64}(undef, m)
#     pi = CuArray{Float64}(undef, m)

#     depind = Int[]
#     sizehint!(depind, m)

#     for i in 1:m
#         # initialize compute components
#         fill!(e_i, 0)
#         e_i[i] = 1

#         # compute pi
#         pi .= A2'[basis,:] \ e_i

#         is_dependent = true
#         for j = 1:n
#             # @show A2[:,j]
#             if !in(j, basis)
#                 view_A2_j = @view A2[:,j]
#                 if abs(view_A2_j' * pi) > 1e-7
#                     basis[i] = j
#                     # println("Pivot: row $i column $j")
#                     is_dependent = false
#                     break
#                 end
#             end
#         end
#         if is_dependent
#             # println("Row $i is dependent.")
#             push!(depind, i)
#         end
#     end
#     # @show basis[basis .> n]
#     # @show basis[basis .> n] .- n

#     if length(depind) > 0
#         deprows = [i for i=1:m if !in(i,depind)]
#         lp.A = lp.A[deprows,:]
#         lp.bl = lp.bl[deprows]
#         lp.bu = lp.bu[deprows]
#         lp.rd = lp.rd[deprows]
#         nrows(lp) -= length(depind)
#         println("Presolve reduced $(length(depind)) rows.")
#     end
# end

"""
Scaling algorithm from Implementation of the Simplex Method, Ping-Qi Pan (2014)
"""
function scaling!(lp::MatOI.AbstractLPForm)
    maxrounds = 50
    round = 1
    while round <= maxrounds
        aratio = matrix_coefficient_ratio(lp.A)
        # @show aratio
        scaling_rows!(lp)

        scaling_columns!(lp)
        sratio = matrix_coefficient_ratio(lp.A)
        # @show sratio, aratio
        if sratio >= 0.9 * aratio
            break
        end
        round += 1
    end
    normalizing_rows!(lp)
    normalizing_columns!(lp)
end

function normalizing_rows!(lp::MatOI.AbstractLPForm{T}) where T
    m, n = size(lp.A)
    rd = zeros(T, m)
    for i = 1:m
        rd[i] = -Inf
        for j = 1:n
            v = abs(lp.A[i,j])
            if v > 0
                rd[i] = rd[i] < v ? v : rd[i]
            end
        end
    end
    lp.A ./= rd
    lp.c_lb ./= rd
    lp.c_ub ./= rd
end

function scaling_rows!(lp::MatOI.AbstractLPForm{T}) where T
    m, n = size(lp.A)
    rd = zeros(T, m)
    for i = 1:m
        mina = Inf; maxa = -Inf
        for j = 1:n
            v = abs(lp.A[i,j])
            if v > 0
                mina = mina > v ? v : mina
                maxa = maxa < v ? v : maxa
            end
        end
        rd[i] = sqrt(mina*maxa)
    end
    lp.A ./= rd
    lp.c_lb ./= rd
    lp.c_ub ./= rd
    println("Scaling rows: max|aij| $(maximum(rd)) min|aij| $(minimum(rd))")
end

function normalizing_columns!(lp::MatOI.AbstractLPForm{T}) where T
    m, n = size(lp.A)
    cd = zeros(T, n)
    vals = nonzeros(lp.A)
    for j = 1:n
        cd[j] = -Inf
        for i in nzrange(lp.A, j)
            v = abs(vals[i])
            if v > 0
                cd[j] = cd[j] < v ? v : cd[j]
            end
        end
    end
    lp.A ./= cd'
    lp.c ./= cd
    lp.v_lb .*= cd
    lp.v_ub .*= cd
end

# function normalizing_columns!(lp::StandardLpData{CuArray})
#     m, n = size(lp.A)
#     for j = 1:n
#         lp.cd[j] = -Inf
#         for i in 1:m
#             v = abs(lp.A[i,j])
#             if v > 0
#                 lp.cd[j] = lp.cd[j] < v ? v : lp.cd[j]
#             end
#         end
#     end
#     lp.A ./= lp.cd'
#     lp.c ./= lp.cd
#     lp.xl .*= lp.cd
#     lp.xu .*= lp.cd
# end

function scaling_columns!(lp::MatOI.AbstractLPForm{T}) where T
    m, n = size(lp.A)
    cd = zeros(T, n)
    vals = nonzeros(lp.A)
    for j = 1:n
        mina = Inf; maxa = -Inf
        for i in nzrange(lp.A, j)
            v = abs(vals[i])
            if v > 0
                mina = mina > v ? v : mina
                maxa = maxa < v ? v : maxa
            end
        end
        cd[j] = sqrt(mina*maxa)
    end
    lp.A ./= cd'
    lp.c ./= cd
    lp.v_lb .*= cd
    lp.v_ub .*= cd
    println("Scaling columns: max|aij| $(maximum(cd)) min|aij| $(minimum(cd))")
end

# function scaling_columns!(lp::StandardLpData{CuArray})
#     m, n = size(lp.A)
#     for j = 1:n
#         mina = Inf; maxa = -Inf
#         for i in 1:m
#             v = abs(lp.A[i,j])
#             if v > 0
#                 mina = mina > v ? v : mina
#                 maxa = maxa < v ? v : maxa
#             end
#         end
#         lp.cd[j] = sqrt(mina*maxa)
#     end
#     lp.A ./= lp.cd'
#     lp.c ./= lp.cd
#     lp.xl .*= lp.cd
#     lp.xu .*= lp.cd
# end

"""
Compute the ratio of the constraint matrix values
"""

function matrix_coefficient_ratio(A::Union{Matrix,CuMatrix})::Float64
    ratio = -Inf
    m, n = size(A)
    for j = 1:n
        mina = Inf; maxa = -Inf
        for i in 1:m
            v = abs(A[i,j])
            if v > 0
                mina = mina > v ? v : mina
                maxa = maxa < v ? v : maxa
            end
        end
        ratio = ratio > maxa / mina ? ratio : maxa / mina
    end
    return ratio
end

function matrix_coefficient_ratio(A::SparseMatrixCSC)::Float64
    ratio = -Inf
    m, n = size(A)
    vals = nonzeros(A)
    for j = 1:n
        mina = Inf; maxa = -Inf
        for i in nzrange(A, j)
            v = abs(vals[i])
            if v > 0
                mina = mina > v ? v : mina
                maxa = maxa < v ? v : maxa
            end
        end
        ratio = ratio > maxa / mina ? ratio : maxa / mina
    end
    return ratio
end


text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot
using ProgressMeter, Statistics
using JLD2

include("koch.jl")

function run_exs_takens(k, F, N, eps_step)

    name = savename("takens_epsilonmax_regularized", @dict(k, N,F,eps_step),"jld2")
    isfile(datadir("takens", name)) && return

    systems = ["Koch", "Henon", "Lorenz", "L96_D5","L96_D6","L96_D7","L96_D8","L96_D9"]
    systs = [pointskoch(points, k), Systems.henon(), Systems.lorenz(),Systems.lorenz96(5, F=F),Systems.lorenz96(6, F=F),Systems.lorenz96(7, F=F),Systems.lorenz96(8, F=F),Systems.lorenz96(9, F=F)]

    Delta_T = Dict()

    for (system, ds) in zip(systems, systs)

        println("System: $system")

        X = nothing

        if system == "Koch"
            points = [[0.0; 0.0], [1.0; 0.0]]
            koch = pointskoch(points, k)
            # making sure we have the same amount of points for the Koch curve
            pt = rand(1:length(koch), 40000)
            X = regularize(Dataset(koch[pt]))
        else
            dt = system == "Henon" ? 1 : 0.1
            T = N * dt
            ts = trajectory(ds, T; dt = dt, Ttr = 100)
            X = regularize(ts)
        end

        εs = 0.2:eps_step:1
        α, αmin, αmax = zeros(length(εs)), zeros(length(εs)), zeros(length(εs))

        @showprogress for (i,ε) in enumerate(εs)
            try
                α[i], αmax[i], αmin[i] = takens_best_estimate(X, ε)
            catch
                # An error above means that the estimator diverged
                α[i], αmax[i], αmin[i] = NaN, NaN, NaN
            end
        end
        result = @strdict α αmin αmax

        Delta_T[system] = result
    end
    @tagsave(datadir("takens", name), Delta_T)
end

function plot_ex_takens(F, N, eps_step,k)

    name = savename("takens_epsilonmax_regularized", @dict(N,F,eps_step,k),"jld2")
    plotname = savename("takens_epsilonmax_regularized_all", @dict(N,F, eps_step,k),"png")

    systems = ["Koch", "Henon", "L96_D5","L96_D7","L96_D9"]
    labels = Dict("Koch"=> "Koch", "Henon"=> "Hénon", "L96_D5"=>"Lorenz-96, \$D = 5\$","L96_D7"=>"Lorenz-96, \$D = 7\$","L96_D9"=>"Lorenz-96, \$D = 9\$")

    @load datadir("takens", name) Koch Henon L96_D5 L96_D7 L96_D9

    system_results = [Koch, Henon, L96_D5, L96_D7, L96_D9]

    @load datadir("lorenz96_Dmax=13_F=$(Int(F))_T=100000_ics=10.jld2") K

    K2 = [log(4)/log(3), 1.26] # analytial solution for Koch, Δ^(L) for Hénon

    Ks = [K2...,mean.(K[[1,3,5]])...]

    εs = 0.2:eps_step:1

    fig, axs = plt.subplots(2, sharex=true,  gridspec_kw=Dict("hspace"=> .1),constrained_layout=false)

    mylinestyles = [":", "--", "-.",":", "--"]

    for (k,(system,system_result)) in enumerate(zip(systems, system_results))

        ax = (system == "Henon" || system == "Koch") ? axs[2] : axs[1]

        ax.fill_between(εs, system_result["αmin"], system_result["αmax"], color = COLORS[k], alpha = .5)
        ax.plot(εs,system_result["α"], color = COLORS[k],linestyle = mylinestyles[k],label = labels[system])

        ax.plot(εs, ones(length(εs))*Ks[k], color = COLORS[k], linestyle = "-", lw = 3)
        ax.set_xlim(.2,1)
    end

    handles, labels = axs[1].get_legend_handles_labels()
    handles2, labels2 = axs[2].get_legend_handles_labels()
    xlabel(L"$\varepsilon_{\mathrm{max}}$")
    fig.text(0.04, 0.55, L"$\Delta^{(T)}$", va="center", rotation="vertical", fontsize=30)
    fig.tight_layout()
    fig.subplots_adjust(bottom=0.13, top = .84)
    fig.legend([handles...,handles2...],[labels...,labels2...], loc = "upper center", ncol = 3)

    tight_layout()
    wsave(plotsdir("takens", plotname), fig)
    close("all")
end

N = 40000
F = 8.
k = 8
run_exs_takens(k,F,N,.01)
plot_ex_takens(F, N,.01,k)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
#=
Protected Accessor Function for User.
=#
# gets the sparse matrix
function getadj(A::Assoc)
    return copy(A.A)
end

# getcol, getrow, getval return range objects if possible
# if not, they return copies of the col, rw, val parts of Assoc's struct
function getcol(A::Assoc)
    if isempty(A)
        return Union{AbstractString, Number}[]
    else
        if isa(A.col[1], AbstractString) || ~(size(A.A,2) == A.col[end])
            return copy(A.col)
        end
        return 1:size(A.A,2)
    end
end

function getrow(A::Assoc)
    if isempty(A)
        return Union{AbstractString, Number}[]
    else
        if isa(A.row[1], AbstractString) || ~(size(A.A,1) == A.row[end])
            return copy(A.row)
        end
        return 1:size(A.A,1)
    end
end

function getval(A::Assoc)
    if isempty(A)
        return Union{AbstractString, Number}[]
    else
        if isa(A.val[1], AbstractString)
            return copy(A.val)
        end
        return sort(unique(nonzeros(A.A)))
    end
end

#=
find : get the triplet of the input Assoc in three array.  Similar to findnz for sparse.
=#
using SparseArrays

function find(A::Assoc)

    if isa(A.A,LinearAlgebra.Adjoint) || isa(A.A, LinearAlgebra.Transpose)
        row, col, val = findnz(SparseMatrixCSC(A.A))
    else
        row, col, val = findnz(A.A)
    end

    val = Array(val)
    n = nnz(A)
    #map if the Associative array isn't numerical
    if ~isempty(A)
        if isa(A.row[1],AbstractString)
            row = [ A.row[row[i]] for i in 1:n]
        end

        if isa(A.col[1],AbstractString)
            col = [ A.col[col[i]] for i in 1:n]
        end

        if isa(A.val[1],AbstractString) || A.val != [1.0]
            val = [ A.val[val[i]] for i in 1:n]
        end
    end

    return row, col, val
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


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

JLLWrappers.@generate_main_file_header("HarfBuzz_ICU")
JLLWrappers.@generate_main_file("HarfBuzz_ICU", UUID("655565e8-fb53-5cb3-b0cd-aec1ca0647ea"))
end  # module HarfBuzz_ICU_jll