text:
# Examples from https://github.com/frewsxcv/python-geojson/tree/master/tests

a = """{
    "properties": {
        "Ã": "Ã"
    },
    "type": "Feature",
    "geometry": null,
    "crs": {
        "type": "link",
        "properties": {
            "href": "data.crs",
            "type": "ogcwkt"
        }
    }
}"""

b = """{"type": "Feature",
    "geometry": {"coordinates": [[-155.52, 19.61], [-156.22, 20.74], [-157.97, 21.46]], "type": "MultiPoint"},
    "id": 1,
    "properties": {"type": "é"}
}"""

c = """{"type": "Feature",
    "geometry": null,
    "id": 1,
    "properties": {"type": "meow"},
    "crs": {"properties": {"name": "urn:ogc:def:crs:EPSG::3785"},
            "type": "name"}
}"""

d = """{
    "type": "Feature",
    "id": "1",
    "bbox": [-180.0, -90.0, 180.0, 90.0],
    "geometry": {"type": "MultiLineString", "coordinates": [[[3.75, 9.25], [-130.95, 1.52]], [[23.15, -34.25], [-1.35, -4.65], [3.45, 77.95]]]},
    "properties": {"title": "Dict 1"}
}"""

e = """{"geometry": {"coordinates": [53, -4],
                     "type": "Point",
                     "crs": {
                         "type": "link",
                         "properties": {
                           "href": "http://example.com/crs/42",
                           "type": "proj4"
                           }
                     }},
    "id": "1",
    "properties": {"link": "http://example.org/features/1",
                   "summary": "The first feature",
                   "title": "Feature 1"},
    "type": "Feature"
}"""

f = """{"geometry": null, "id": 12, "properties": {"foo": "bar"}, "type": "Feature"}"""

# Examples from https://github.com/Esri/geojson-utils/blob/master/tests/geojson.js

multipolygon = """{
    "type": "FeatureCollection",
    "features": [{
        "type": "Feature",
        "bbox": [100.0, 0.0, 105.0, 1.0],
        "properties": {
            "cartodb_id": 46,
            "addr1": "18150 E. Pathfinder Rd.",
            "addr2": "Rowland Heights",
            "park": "Pathfinder Park"
        },
        "geometry": {
            "type":"MultiPolygon",
            "coordinates": [[[ [-117.913883,33.96657], [-117.907767,33.967747], [-117.912919,33.96445], [-117.913883,33.96657] ]]]
        }
    }]
}"""

realmultipolygon = """{
        "type":"FeatureCollection",
        "features":[
        {
            "type":"Feature",
            "properties":{
                "cartodb_id":46,
                "addr1":"18150 E. Pathfinder Rd.",
                "addr2":"Rowland Heights",
                "park":"Pathfinder Park"
            },
            "geometry":{
                "type":"MultiPolygon",
                "coordinates": [
                    [[[102.0, 2.0], [103.0, 2.0], [103.0, 3.0], [102.0, 3.0], [102.0, 2.0]]],
                [[[100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0]],
                [[100.2, 0.2], [100.8, 0.2], [100.8, 0.8], [100.2, 0.8], [100.2, 0.2]]]
                    ]
            }
        }
        ]
    }"""

polyline = """{
        "next_feature" : "1",
        "type" : "FeatureCollection",
        "start" : 0,
        "features" : [{
            "type" : "Feature",
            "id" : "a73ws67n775q",
            "geometry" : {
                "type" : "LineString",
                "coordinates" : [[-89, 43], [-88, 44], [-88, 45]]
            },
            "properties" : {
                "InLine_FID" : 0,
                "SimLnFLag" : 0
            }
        }],
        "sort" : null,
        "page" : 0,
        "count" : 2073,
        "limit" : 1
    }"""

point = """{
        "next_feature" : "1",
        "type" : "FeatureCollection",
        "start" : 0,
        "features" : [{
            "type" : "Feature",
            "id" : "a7vs0i9rnyyx",
            "geometry" : {
                "type" : "Point",
                "coordinates" : [-89, 44]
            },
            "properties" : {
                "fax" : "305-571-8347",
                "phone" : "305-571-8345"
            }
        }],
        "sort" : null,
        "page" : 0,
        "count" : 236,
        "limit" : 1
    }"""

pointnull = """{
        "next_feature" : "1",
        "type" : "FeatureCollection",
        "start" : 0,
        "features" : [{
            "type" : "Feature",
            "id" : "a7vs0i9rnyyx",
            "geometry" : null,
            "properties" : {
                "fax" : "305-571-8347",
                "phone" : "305-571-8345"
            }
        }],
        "sort" : null,
        "page" : 0,
        "count" : 236,
        "limit" : 1
    }"""

poly = """{
        "next_feature" : "1",
        "type" : "FeatureCollection",
        "start" : 0,
        "features" : [{
            "type" : "Feature",
            "id" : "a7ws7wldxold",
            "geometry" : {
                "type" : "Polygon",
                "coordinates" : [[[-89, 42], [-89, 50], [-80, 50], [-80, 42], [-89, 42]]]
            },
            "properties" : {
                "DIST_NUM" : 7.0,
                "LOCATION" : "Bustleton Ave. & Bowler St",
                "PHONE" : "686-3070",
                "DIST_NUMC" : "07",
                "DIV_CODE" : "NEPD",
                "AREA_SQMI" : 12.41643
            }
        }],
        "sort" : null,
        "page" : 0,
        "count" : 25,
        "limit" : 1
    }"""

polyhole = """{
        "next_feature" : "1",
        "type" : "FeatureCollection",
        "start" : 0,
        "features" : [{
            "type" : "Feature",
            "id" : "a7ws7wldxold",
            "geometry" : {
                "type" : "Polygon",
                "coordinates" : [[[-89, 42], [-89, 50], [-80, 50], [-80, 42], [-89, 42]], [[-87, 44], [-82, 44], [-82, 48], [-87, 48], [-87, 44]]]
            },
            "properties" : {
                "DIST_NUM" : 7.0,
                "LOCATION" : "Bustleton Ave. & Bowler St",
                "PHONE" : "686-3070",
                "DIST_NUMC" : "07",
                "DIV_CODE" : "NEPD",
                "AREA_SQMI" : 12.41643
            }
        }],
        "sort" : null,
        "page" : 0,
        "count" : 25,
        "limit" : 1
    }"""

collection = """{
        "next_feature" : "1",
        "type" : "FeatureCollection",
        "start" : 0,
        "features" : [{
            "type" : "Feature",
            "id" : "a7xlmuwyjioy",
            "geometry" : {
                "type" : "GeometryCollection",
                "geometries" : [{
                    "type" : "Polygon",
                    "coordinates" : [[[-95, 43], [-95, 50], [-90, 50], [-91, 42], [-95, 43]]]
                }, {
                    "type" : "Polygon",
                    "coordinates" : [[[-89, 42], [-89, 50], [-80, 50], [-80, 42], [-89, 42]]]
                }, {
                    "type" : "Point",
                    "coordinates" : [-94, 46]
                }]
            },
            "properties" : {
                "STATE_ABBR" : "ZZ",
                "STATE_NAME" : "Top"
            }
        }],
        "sort" : null,
        "page" : 0,
        "count" : 3,
        "limit" : 1
    }"""

# Example from osmbuildings.org/examples/GeoJSON.php

osm_buildings = """{
  "type": "FeatureCollection",
  "features": [{
    "type": "Feature",
    "geometry": {
      "type": "Polygon",
      "coordinates": [
        [
          [13.42634, 52.49533],
          [13.42660, 52.49524],
          [13.42619, 52.49483],
          [13.42583, 52.49495],
          [13.42590, 52.49501],
          [13.42611, 52.49494],
          [13.42640, 52.49525],
          [13.42630, 52.49529],
          [13.42634, 52.49533]
        ]
      ]
    },
    "properties": {
      "color": "rgb(255,200,150)",
      "height": 150
    }
  }, {
    "type": "Feature",
    "geometry": {
      "type": "Polygon",
      "coordinates": [
        [
          [13.42706, 52.49535],
          [13.42745, 52.49520],
          [13.42745, 52.49520],
          [13.42741, 52.49516],
          [13.42717, 52.49525],
          [13.42692, 52.49501],
          [13.42714, 52.49494],
          [13.42686, 52.49466],
          [13.42650, 52.49478],
          [13.42657, 52.49486],
          [13.42678, 52.49480],
          [13.42694, 52.49496],
          [13.42675, 52.49503],
          [13.42706, 52.49535]
        ]
      ]
    },
    "properties": {
      "color": "rgb(180,240,180)",
      "height": 130
    }
  }, {
    "type": "Feature",
    "geometry": {
      "type": "MultiPolygon",
      "coordinates": [
        [
          [
            [13.42746, 52.49440],
            [13.42794, 52.49494],
            [13.42799, 52.49492],
            [13.42755, 52.49442],
            [13.42798, 52.49428],
            [13.42846, 52.49480],
            [13.42851, 52.49478],
            [13.42800, 52.49422],
            [13.42746, 52.49440]
          ]
        ],
        [
          [
            [13.42803, 52.49497],
            [13.42800, 52.49493],
            [13.42844, 52.49479],
            [13.42847, 52.49483],
            [13.42803, 52.49497]
          ]
        ]
      ]
    },
    "properties": {
      "color": "rgb(200,200,250)",
      "height": 120
    }
  }, {
    "type": "Feature",
    "geometry": {
      "type": "Polygon",
      "coordinates": [
        [
          [13.42857, 52.49480],
          [13.42918, 52.49465],
          [13.42867, 52.49412],
          [13.42850, 52.49419],
          [13.42896, 52.49465],
          [13.42882, 52.49469],
          [13.42837, 52.49423],
          [13.42821, 52.49428],
          [13.42863, 52.49473],
          [13.42853, 52.49476],
          [13.42857, 52.49480]

        ]
      ]
    },
    "properties": {
      "color": "rgb(150,180,210)",
      "height": 140
    }
  }]
}"""


text:
# [MLJ.jl]: https://github.com/alan-turing-institute/MLJ.jl
# [RDatasets.jl]: https://github.com/JuliaStats/RDatasets.jl
# [MLJModels.jl]: https://github.com/alan-turing-institute/MLJModels.jl
# [DecisionTree.jl]: https://github.com/bensadeghi/DecisionTree.jl
# [NearestNeighbors.jl]: https://github.com/KristofferC/NearestNeighbors.jl
# [GLM.jl]: https://github.com/JuliaStats/GLM.jl
# [ScikitLearn.jl]: https://github.com/cstjean/ScikitLearn.jl
# ## Data and its interpretation
#
# ### Machine type and scientific type

using RDatasets, MLJ
MLJ.color_off() # hide
iris = dataset("datasets", "iris")

first(iris, 3) |> pretty

# Observe that below each column name there are two _types_ given: the first one is the _machine type_ and the second one is the _scientific type_.
#
# * **machine type**: is the Julia type the data is currently encoded as, for instance `Float64`,
# * **scientific type**: is a type corresponding to how the data should be _interpreted_, for instance `Multiclass{3}`.
#
# If you want to specify a different scientific type than the one inferred, you can do so by using the function `coerce` along with pairs of column names and scientific types:

iris2 = coerce(iris, :PetalWidth => OrderedFactor)
first(iris2[:, [:PetalLength, :PetalWidth]], 1) |> pretty

# ### Unpacking data
#
# The function `unpack` helps specify the target and the input for a regression or classification task

y, X = unpack(iris, ==(:Species), colname -> true)
first(X, 1) |> pretty

# The two arguments after the dataframes should be understood as _functions_ over column names specifying the target and the input data respectively.
# Let's look in more details at what we used here:
#
# * `==(:Species)` is a shorthand to specify that the target should be the column with name equal to `:Species`,
# * `colname -> true` indicates that every other column is to be taken as input
#
# Let's try another one:

y, X = unpack(iris, ==(:Species), !=(:PetalLength))
first(X, 1) |> pretty

# You can also use the shorthand `@load_iris` for such common examples:

X, y = @load_iris;

# ## Choosing a model
#
# ### Model search
#
# In MLJ, a _model_ is a struct storing the _hyperparameters_ of the learning algorithm indicated by the struct name (and only that).
#
# A number of models are available in MLJ, usually thanks to external packages interfacing with MLJ (see also `[MLJModels.jl]`).
# In order to see which ones are appropriate for the data you have and its scientific interpretation, you can use the function `models` along with the function `matching`; let us look specifically at models which support a probabilistic output:

for m in models(matching(X, y))
    if m.prediction_type == :probabilistic
        println(rpad(m.name, 30), "($(m.package_name))")
    end
end

# ### Loading a model
#
# Most models are implemented outside of the MLJ ecosystem; you therefore have to _load models_ using the `@load` command.
#
# **Note**: you _must_ have the package from which the model is loaded available in your environment (in this case `[DecisionTree.jl]`) otherwise MLJ will not be able to load the model.
#
# For instance, let's say you want to fit a K-Nearest Neighbours classifier:

knc = @load KNeighborsClassifier

# In some cases, there may be several packages offering the same model, for instance `LinearRegressor` is offered by both `[GLM.jl]` and `[ScikitLearn.jl]` so you will need to specify the package you would like to use by adding `pkg="ThePackage"` in the load command:

linreg = @load LinearRegressor pkg="GLM"

ext: jl

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
function rotatelattice(A, B, C)

Anorm = norm(A);
Bnorm = norm(B);
Cnorm = norm(C);
Ahat = A/Anorm;

ax = Anorm;
bx = dot(B,Ahat);
by = sqrt(Bnorm^2 - bx^2); #norm(cross(Ahat,B));
cx = dot(C,Ahat);
cy = (dot(B, C) - bx*cx)/by;
cz = sqrt(Cnorm^2 - cx^2 - cy^2);

a = [ax 0 0];
b = [bx by 0];
c = [cx cy cz];

return a, b, c

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# * BitFloats

module BitFloats

export Float80,  Inf80,  NaN80,
       Float128, Inf128, NaN128

import Base: !=, *, +, -, /, <, <=, ==, ^, abs, bswap, cos, decompose, eps, exp, exp2,
             exponent, exponent_half, exponent_mask, exponent_one, floatmax, floatmin,
             isequal, isless, issubnormal, ldexp, log, log10, log2, nextfloat, precision,
             promote_rule, reinterpret, rem, round, show, sign_mask, significand,
             significand_mask, sin, sqrt, trunc, typemax, typemin, uinttype, unsafe_trunc

import Base.Math: exponent_bits, exponent_raw_max, significand_bits

using Base: bswap_int, llvmcall, uniontypes

using Base.GMP: BITS_PER_LIMB, Limb

using Core: BuiltinInts

import BitIntegers

import Random: rand
using  Random: AbstractRNG, CloseOpen01, CloseOpen12, SamplerTrivial

# * definitions

primitive type Float80  <: AbstractFloat 80  end
primitive type Float128 <: AbstractFloat 128 end

const WBF = Union{Float80,Float128}

BitIntegers.@define_integers 80

uinttype(::Type{Float80}) = UInt80
uinttype(::Type{Float128}) = UInt128

const llvmvars = ((Float80, "x86_fp80", "i80", "f80"), (Float128, "fp128", "i128", "f128"))

# * traits

sign_mask(::Type{Float80})            = 0x8000_0000_0000_0000_0000 % UInt80
exponent_mask(::Type{Float80})        = 0x7fff_0000_0000_0000_0000 % UInt80
exponent_one(::Type{Float80})         = 0x3fff_8000_0000_0000_0000 % UInt80
exponent_half(::Type{Float80})        = 0x3ffe_8000_0000_0000_0000 % UInt80
significand_mask(::Type{Float80})     = 0x0000_ffff_ffff_ffff_ffff % UInt80
explicit_bit(::Type{Float80}=Float80) = 0x0000_8000_0000_0000_0000 % UInt80
# non-implicit most significand bit of significand is actually stored for Float80

sign_mask(::Type{Float128})           = 0x8000_0000_0000_0000_0000_0000_0000_0000
exponent_mask(::Type{Float128})       = 0x7fff_0000_0000_0000_0000_0000_0000_0000
exponent_one(::Type{Float128})        = 0x3fff_0000_0000_0000_0000_0000_0000_0000
exponent_half(::Type{Float128})       = 0x3ffe_0000_0000_0000_0000_0000_0000_0000
significand_mask(::Type{Float128})    = 0x0000_ffff_ffff_ffff_ffff_ffff_ffff_ffff

significand_bits(::Type{T}) where {T<:WBF} = trailing_ones(significand_mask(T))
exponent_bits(   ::Type{T}) where {T<:WBF} = sizeof(T)*8 - significand_bits(T) - 1
exponent_bias(   ::Type{T}) where {T<:WBF} = Int(exponent_one(T) >> significand_bits(T))
# exponent_bias is 16383 for both types
exponent_raw_max(::Type{T}) where {T<:WBF} = Int(exponent_mask(T) >> significand_bits(T))

eps(     ::Type{Float80})  = reinterpret(Float80,  0x3fc0_8000_0000_0000_0000 % UInt80)
floatmin(::Type{Float80})  = reinterpret(Float80,  0x0001_8000_0000_0000_0000 % UInt80)
floatmax(::Type{Float80})  = reinterpret(Float80,  0x7ffe_ffff_ffff_ffff_ffff % UInt80)
typemin( ::Type{Float80})  = -Inf80
typemax( ::Type{Float80})  =  Inf80

eps(     ::Type{Float128}) = reinterpret(Float128, 0x3f8f_0000_0000_0000_0000_0000_0000_0000)
floatmin(::Type{Float128}) = reinterpret(Float128, 0x0001_0000_0000_0000_0000_0000_0000_0000)
floatmax(::Type{Float128}) = reinterpret(Float128, 0x7ffe_ffff_ffff_ffff_ffff_ffff_ffff_ffff)
typemin( ::Type{Float128}) = -Inf128
typemax( ::Type{Float128}) =  Inf128

const Inf80 = reinterpret(Float80,   0x7fff_8000_0000_0000_0000 % UInt80)
const NaN80 = reinterpret(Float80,   0x7fff_c000_0000_0000_0000 % UInt80)

const Inf128 = reinterpret(Float128, 0x7fff_0000_0000_0000_0000_0000_0000_0000)
const NaN128 = reinterpret(Float128, 0x7fff_8000_0000_0000_0000_0000_0000_0000)

precision(::Type{Float80})  = 64
precision(::Type{Float128}) = 113

# * float functions

function exponent(x::T) where T<:WBF
    @noinline throw1(x) = throw(DomainError(x, "Cannot be NaN or Inf."))
    @noinline throw2(x) = throw(DomainError(x, "Cannot be subnormal converted to 0."))
    xs = reinterpret(Unsigned, x) & ~sign_mask(T)
    xs >= exponent_mask(T) && throw1(x)
    k = Int(xs >> significand_bits(T))
    if k == 0 # x is subnormal
        xs == 0 && throw2(x)
        m = leading_zeros(xs) - exponent_bits(T)
        if x isa Float80
            m -= 1 # non-implicit bit
        end
        k = 1 - m
    end
    return k - exponent_bias(T)
end

function significand(x::T) where T<:WBF
    xu = reinterpret(Unsigned, x)
    xs = xu & ~sign_mask(T)
    xs >= exponent_mask(T) && return x # NaN or Inf
    if xs <= significand_mask(T) # x is subnormal
        xs == 0 && return x # +-0
        m = unsigned(leading_zeros(xs) - exponent_bits(T)) - (x isa Float80)
        xs <<= m
        xu = xs | (xu & sign_mask(T))
        if x isa Float80
            xu |= explicit_bit() # set non-implicit bit to 1
        end
    end
    xu = (xu & ~exponent_mask(T)) | exponent_one(T)
    return reinterpret(T, xu)
end

function squeezeimplicit(u::UInt80)
    s = ((u & significand_mask(Float80)) << 1) & significand_mask(Float80)
    ((u & ~significand_mask(Float80)) | s) >> 1
end

function unsqueezeimplicit(u::UInt80)
    e = (u & ~(significand_mask(Float80) >> 1)) << 1
    u = e | (u & (significand_mask(Float80) >> 1))
    iszero(exponent_mask(Float80) & u) ? # zero or subnormal
        u :
        u | explicit_bit()
end

squeezeimplicit(f::UInt128) = f
unsqueezeimplicit(f::UInt128) = f

function nextfloat(f::WBF, d::Integer)
    F = typeof(f)
    fumax = squeezeimplicit(reinterpret(Unsigned, F(Inf)))
    U = typeof(fumax)

    isnan(f) && return f
    fi = reinterpret(Signed, f)
    fneg = fi < 0
    fu = squeezeimplicit(unsigned(fi & typemax(fi)))

    dneg = d < 0
    da = Base.uabs(d)
    if da > typemax(U)
        fneg = dneg
        fu = fumax
    else
        du = da % U
        if fneg ⊻ dneg
            if du > fu
                fu = min(fumax, du - fu)
                fneg = !fneg
            else
                fu = fu - du
            end
        else
            if fumax - fu < du
                fu = fumax
            else
                fu = fu + du
            end
        end
    end
    fu = unsqueezeimplicit(fu)
    if fneg
        fu |= sign_mask(F)
    end
    reinterpret(F, fu)
end

# for Float80, left-most (non-implicit) bit of significand should be 0
function issubnormal(x::T) where {T<:WBF}
    y = reinterpret(Unsigned, x)
    (y & exponent_mask(T) == 0) & (y & significand_mask(T) != 0)
end

function ldexp(x::T, e::Integer) where T<:WBF
    xu = reinterpret(Unsigned, x)
    xs = xu & ~sign_mask(T)
    xs >= exponent_mask(T) && return x # NaN or Inf
    k = Int(xs >> significand_bits(T))
    if k == 0 # x is subnormal
        xs == 0 && return x # +-0
        m = leading_zeros(xs) - exponent_bits(T) - (x isa Float80)
        ys = xs << unsigned(m)
        xu = ys | (xu & sign_mask(T))
        k = 1 - m
        # underflow, otherwise may have integer underflow in the following n + k
        e < -50000 && return flipsign(T(0.0), x)
    end
    # For cases where e of an Integer larger than Int make sure we properly
    # overflow/underflow; this is optimized away otherwise.
    if e > typemax(Int)
        return flipsign(T(Inf), x)
    elseif e < typemin(Int)
        return flipsign(T(0.0), x)
    end
    n = e % Int
    k += n
    # overflow, if k is larger than maximum possible exponent
    if k >= exponent_raw_max(T)
        return flipsign(T(Inf), x)
    end
    if k > 0 # normal case
        xu = (xu & ~exponent_mask(T)) | (rem(k, uinttype(T)) << significand_bits(T))
        if x isa Float80
            xu |= explicit_bit() # was
        end
        return reinterpret(T, xu)
    else # subnormal case
        expo = significand_bits(T) - (x isa Float80)
        if k <= -expo # underflow
            # overflow, for the case of integer overflow in n + k
            e > 50000 && return flipsign(T(Inf), x)
            return flipsign(T(0.0), x)
        end
        k += expo # k > 0
        z = T(2.0)^-expo
        xu = (xu & ~exponent_mask(T)) | (rem(k, uinttype(T)) << significand_bits(T))
        x isa Float80 && (xu |= explicit_bit())
        return z*reinterpret(T, xu)
    end
end

# * conversions

# ** signedness

reinterpret(::Type{Unsigned}, x::Float80)  = reinterpret(UInt80, x)
reinterpret(::Type{Unsigned}, x::Float128) = reinterpret(UInt128, x)
reinterpret(::Type{Signed},   x::Float80)  = reinterpret(Int80, x)
reinterpret(::Type{Signed},   x::Float128) = reinterpret(Int128, x)

Signed(  x::WBF) = Int(x)
Unsigned(x::WBF) = UInt(x)

trunc(::Type{Signed},   x::WBF) = trunc(Int, x)
trunc(::Type{Integer},  x::WBF) = trunc(Int, x)
trunc(::Type{Unsigned}, x::WBF) = trunc(UInt, x)

# ** from ints

# unsafe_trunc & F(::T)
for (F, f, i) = llvmvars
    for T = uniontypes(BuiltinInts)
        s = 8*sizeof(T)

        # T -> F
        itofp = T <: Signed ? :sitofp : :uitofp
        @eval begin
            (::Type{$F})(x::$T) = llvmcall(
                $"""
                %y = $itofp i$s %0 to $f
                %mi = bitcast $f %y to $i
                ret $i %mi
                """,
                $F, Tuple{$T}, x)

            promote_rule(::Type{$F}, ::Type{$T}) = $F
        end
        T === Bool && continue

        # F -> T
        fptoi = T <: Signed ? :fptosi : :fptoui
        @eval unsafe_trunc(::Type{$T}, x::$F) = llvmcall(
            $"""
            %x = bitcast $i %0 to $f
            %y = $fptoi $f %x to i$s
            ret i$s %y
            """,
            $T, Tuple{$F}, x)
    end
end

# trunc & T(::F)
for Tf in (Float80, Float128)
    for Ti in Base.BitInteger_types
        Ti == Int128 && continue
        @eval begin
            function trunc(::Type{$Ti},x::$Tf)
                # unlike in Base, these expressions can't be computed at compile time,
                # this segfaults otherwise
                if $Tf(typemin($Ti)) - one($Tf) < x < $Tf(typemax($Ti)) + one($Tf)
                    return unsafe_trunc($Ti, x)
                else
                    throw(InexactError(:trunc, $Ti, x))
                end
            end

            function (::Type{$Ti})(x::$Tf)
                if $Tf(typemin($Ti)) <= x < $Tf(typemax($Ti)) + one($Tf) && (round(x, RoundToZero) == x)
                    return unsafe_trunc($Ti, x)
                else
                    throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
                end
            end
        end
    end

    Ti = Int128
    @eval begin
        function trunc(::Type{$Ti},x::$Tf)
            if $Tf(typemin($Ti)) <= x < $Tf(typemax($Ti))
                return unsafe_trunc($Ti,x)
            else
                throw(InexactError(:trunc, $Ti, x))
            end
        end

        function (::Type{$Ti})(x::$Tf)
            if ($Tf(typemin($Ti)) <= x < $Tf(typemax($Ti))) && (round(x, RoundToZero) == x)
                return unsafe_trunc($Ti,x)
            else
                throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
            end
        end
    end
end

# ** from floats

for (F, f, i) = llvmvars
    for (S, s) = ((Float32, :float), (Float64, :double))
        @eval begin
            (::Type{$F})(x::$S) = llvmcall(
                $"""
                %y = fpext $s %0 to $f
                %yi = bitcast $f %y to $i
                ret $i %yi
                """,
                $F, Tuple{$S}, x)

            (::Type{$S})(x::$F) = llvmcall(
                $"""
                %x = bitcast $i %0 to $f
                %y = fptrunc $f %x to $s
                ret $s %y
                """,
                $S, Tuple{$F}, x)

            promote_rule(::Type{$F}, ::Type{$S}) = $F
        end
    end
    @eval begin
        (::Type{$F})(x::Float16) = $F(Float32(x))
        (::Type{Float16})(x::$F) = Float16(Float32(x))
        promote_rule(::Type{$F}, ::Type{Float16}) = $F
    end
end

# ** Float80 <-> Float128

promote_rule(::Type{Float128}, ::Type{Float80}) = Float128

if false # unimplemented by llvm
    @eval begin
        (::Type{Float128})(x::Float80) = llvmcall(
            """
            %x = bitcast i80 %0 to x86_fp80
            %y = fpext x86_fp80 %x to fp128
            %yi = bitcast fp128 %y to i128
            ret i128 %yi
            """,
            Float128, Tuple{Float80}, x)

        (::Type{Float80})(x::Float128) = llvmcall(
            """
            %x = bitcast i128 %0 to fp128
            %y = fptrunc fp128 %x to x86_fp80
            %yi = bitcast x86_fp80 %y to i80
            ret i80 %yi
            """,
            Float80, Tuple{Float128}, x)
    end
else
    function (::Type{Float128})(x::Float80)
        u = reinterpret(Unsigned, x)
        se = ((u & ~significand_mask(Float80)) % UInt128) << 48
        u &= significand_mask(Float80) >> 1 # >> 1 for explicit bit
        reinterpret(Float128, se | ((u % UInt128) << 49))
    end

    function (::Type{Float80})(x::Float128)
        u = reinterpret(Unsigned, x)
        se = ((u & ~significand_mask(Float128)) >> 48) % UInt80
        v = ((u & significand_mask(Float128)) >> 49) % UInt80
        if se & exponent_mask(Float80) != 0
            v |= explicit_bit()
        end
        reinterpret(Float80, se | v)
    end
end

# ** BigFloat

function (::Type{BigFloat})(x::WBF; precision::Integer=Base.MPFR.DEFAULT_PRECISION[])
    if isnan(x)
        BigFloat(NaN, precision=precision)
    elseif isinf(x)
        BigFloat(x < 0 ? -Inf : Inf, precision=precision)
    elseif iszero(x)
        BigFloat(x < 0 ? -0.0 : 0.0, precision=precision)
    else
        z = BigFloat(precision=precision)
        z.exp = exponent(x) + 1
        z.sign = x >= zero(x) ? 1 : -1
        y = significand(x)
        nlimbs = div(precision + BITS_PER_LIMB-1, BITS_PER_LIMB)
        mask = ~(Limb(0)) << mod(BITS_PER_LIMB - precision, BITS_PER_LIMB)

        if x isa Float80
            u = reinterpret(Unsigned, y) % UInt64
            n = nlimbs - div(64, BITS_PER_LIMB) + 1
        else
            u = reinterpret(Unsigned, y) % UInt128
            u = (u << 16) >> 1
            u |= sign_mask(Float128) # explicit 1 for the msb
            n = nlimbs - div(128, BITS_PER_LIMB) + 1
        end

        for i = 1:nlimbs
            l = (u >> ((i-n) * BITS_PER_LIMB)) % Limb
            i == 1 && (l &= mask)
            GC.@preserve z unsafe_store!(z.d, l, i)
        end
        z
    end
end

# alternative for Float128, useful to compare accuracy
function BigFloat_mpfr(x::Float128)
    z = BigFloat()
    ccall((:mpfr_set_float128, :libmpfr), Int32, (Ref{BigFloat}, Float128, Int32), z, x, Base.MPFR.ROUNDING_MODE[])
    if isnan(x) && signbit(x) != signbit(z)
        z.sign = -z.sign
    end
    return z
end

function Float128(x::BigFloat)
    ispos = x.sign == 1
    if isnan(x)
        NaN128
    elseif iszero(x) # should not be necessary (i.e. it works in the else branch too)
       x.sign * zero(Float128)
    else
        expo = x.exp - 1 # exponent bias for mpfr is -1
        # comparing x > floatmax(Float128) allocates a BigFloat, so we look at
        # expo to check for infinity
        if expo > 16383 || isinf(x)
            return x.sign * Inf128
        end
        e = Int(expo + exponent_bias(Float128))
        @assert e <= Int(0x7fff)
        n = ((x.prec-1) >> ( 6 - (BITS_PER_LIMB === 32))) + 1 # should be the number of limbs
        @assert n >= 1
        GC.@preserve x if BITS_PER_LIMB === 64
            u = (unsafe_load(x.d, n) % UInt128) << 64
            if n > 1
                u |= unsafe_load(x.d, n-1)
            end
        else
            u = zero(UInt128)
            shift = 96
            while n > 0 && shift >= 0
                u |= (unsafe_load(x.d, n) % UInt128) << shift
                n -= 1
                shift -= 32
            end
        end
        if e <= 0 # subnormal
            u >>= -e
            e = 0
        else
            u <<= 1 # remove explicit bit
        end
        u >>= exponent_bits(Float128) + 1
        u |= UInt128(e) << significand_bits(Float128)
        u |= sign_mask(Float128) * !ispos
        reinterpret(Float128, u)
    end
end

# doesn't work properly :(
# Float128_mpfr(x::BigFloat) =
#    ccall((:mpfr_get_float128, :libmpfr), Float128, (Ref{BigFloat}, Int32), x, Base.MPFR.ROUNDING_MODE[])

Float80(x::BigFloat) = Float80(Float128(x))

# ** round

for (F, f, i, fn) = llvmvars
    # TODO: can be broken for Float128
    for (mode, llvmfun) = ((:ToZero, :trunc), (:Down, :floor),
                           (:Up, :ceil), (:Nearest, :rint))
        fun = "@llvm.$llvmfun.$fn"
        @eval round(x::$F, r::$(RoundingMode{mode})) = llvmcall(
            ($"""declare $f $fun($f %Val)""",
             $"""
             %x = bitcast $i %0 to $f
             %y = call $f $fun($f %x)
             %z = bitcast $f %y to $i
             ret $i %z
             """), $F, Tuple{$F}, x)
    end
end

# * comparisons

for (F, f, i) = llvmvars
    for (op, fop) = ((:(==), :oeq), (:!=, :une), (:<, :olt), (:<=, :ole))
        @eval $op(x::$F, y::$F) = llvmcall(
            $"""
            %x = bitcast $i %0 to $f
            %y = bitcast $i %1 to $f
            %b = fcmp $fop $f %x, %y
            %c = zext i1 %b to i8
            ret i8 %c
            """,
            Bool, Tuple{$F,$F}, x, y)
    end

    # adapted from fpislt in intrinsics.cpp
    # could as well be written in Julia
    @eval fpislt(x::$F, y::$F) = llvmcall(
        $"""
        ; %xi = %0
        ; %yi = %1
        %xf = bitcast $i %0 to $f
        %yf = bitcast $i %1 to $f
        %c11 = fcmp ord $f %xf, %xf ; !isnan(xf)
        %c12 = fcmp uno $f %yf, %yf ; isnan(yf)
        %c1  = or i1 %c11, %c12
        %c21 = fcmp ord $f %xf, %yf ; !(isnan(xf) | isnan(yf))
        %c2211 = icmp sge $i %0, 0 ; xi >= 0
        %c2212 = icmp slt $i %0, %1 ; xi < yi
        %c221 = and i1 %c2211, %c2212
        %c2221 = icmp slt $i %0, 0 ; xi < 0
        %c2222 = icmp ugt $i %0, %1 ; unsigned(xi) > unsigned(yi)
        %c222 = and i1 %c2221, %c2222
        %c22 = or i1 %c221, %c222
        %c2 = and i1 %c21, %c22
        %c = or i1 %c1, %c2
        %cbool = zext i1 %c to i8
        ret i8 %cbool
        """,
        Bool, Tuple{$F,$F}, x, y)
end

isless(x::T, y::T) where {T<:WBF} = fpislt(x, y)

# implements Base.fpiseq
isequal(x::T, y::T) where {T<:WBF} =
    (isnan(x) & isnan(y)) | (reinterpret(Unsigned, x) === reinterpret(Unsigned, y))

# * arithmetic

for (F, f, i, fn) = llvmvars
    for (op, fop) = ((:*, :fmul), (:/, :fdiv), (:+, :fadd), (:-, :fsub), (:rem, :frem))
        @eval $op(x::$F, y::$F) = llvmcall(
            $"""
            %x = bitcast $i %0 to $f
            %y = bitcast $i %1 to $f
            %m = $fop $f %x, %y
            %mi = bitcast $f %m to $i
            ret $i %mi
            """, $F, Tuple{$F,$F}, x, y)
    end
    for (op, fop) = (:abs => :fabs, :log2 => :log2, :exp2 => :exp2, :sqrt => :sqrt,
                     :sin => :sin, :cos => :cos, :exp => :exp, :log => :log, :log10 => :log10)
        if F === Float128 && op ∈ (:log2, :exp2, :sqrt, :sin, :cos, :exp, :log, :log10)
            @eval $op(x::$F) = $F($op(big(x)))
        else
            @eval $op(x::$F) = llvmcall(
                ($"""declare $f  @llvm.$fop.$fn($f %Val)""",
                 $"""
                 %x = bitcast $i %0 to $f
                 %y = call $f @llvm.$fop.$fn($f %x)
                 %z = bitcast $f %y to $i
                 ret $i %z
                 """), $F, Tuple{$F}, x)
        end
    end
    F == Float128 && continue # broken
    @eval ^(x::$F, y::$F) = llvmcall(
        ($"""declare $f @llvm.pow.$fn($f %Val, $f %Power)""",
         $"""
         %x = bitcast $i %0 to $f
         %y = bitcast $i %1 to $f
         %m = call $f @llvm.pow.$fn($f %x, $f %y)
         %mi = bitcast $f %m to $i
         ret $i %mi
         """), $F, Tuple{$F,$F}, x, y)
end

^(x::Float128, y::Float128) = Float128(big(x)^big(y))

^(x::WBF, n::Integer) = x^(oftype(x, n))
-(x::F) where {F<:WBF} = F(-0.0) - x

# * hashing

# could return UInt64 instead of Int128 for Float80
function decompose(x::WBF)::Tuple{Int128,Int,Int}
    T = typeof(x)
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(Unsigned, x)
    s = (n & significand_mask(T)) % Int128
    e = ((n & exponent_mask(T)) >> significand_bits(T)) % Int
    if T === Float128 # add explicit bit
        s |= Int128(e != 0) << significand_bits(T)
    end
    d = ifelse(signbit(x), -1, 1)
    # for Float80, significand_bits(T) includes the explicit bit
    E = exponent_bias(T) + significand_bits(T) - (T === Float80)
    s, e - E + (e == 0), d
end

# * rand

rand(rng::AbstractRNG, ::SamplerTrivial{CloseOpen12{Float80}}) =
    reinterpret(Float80, rand(rng, UInt64) | exponent_one(Float80) | explicit_bit())

rand(rng::AbstractRNG, ::SamplerTrivial{CloseOpen12{Float128}}) =
    reinterpret(Float128, rand(rng, UInt128) & significand_mask(Float128) | exponent_one(Float128))

rand(rng::AbstractRNG, ::SamplerTrivial{CloseOpen01{T}}) where {T<:WBF} =
    rand(rng, CloseOpen12(T)) - one(T)

# * misc

bswap(x::WBF) = bswap_int(x)

show(io::IO, x::WBF) =
    if isnan(x)
        print(io, "NaN", sizeof(x)*8)
    elseif isinf(x)
        print(io, x > 0 ? "" : '-',  "Inf", sizeof(x)*8)
    else
        show(io, BigFloat(x, precision=precision(x)))
    end

end # module


text:
export SubscribableTrait, SimpleSubscribableTrait, ScheduledSubscribableTrait, InvalidSubscribableTrait
export AbstractSubscribable, Subscribable, ScheduledSubscribable, as_subscribable
export subscribe!, on_subscribe!
export subscribable_extract_type

export InvalidSubscribableTraitUsageError, InconsistentActorWithSubscribableDataTypesError
export MissingOnSubscribeImplementationError, MissingOnScheduledSubscribeImplementationError

import Base: show, showerror
import Base: eltype

"""
Abstract type for all possible subscribable traits

See also: [`SimpleSubscribableTrait`](@ref), [`ScheduledSubscribableTrait`](@ref), [`InvalidSubscribableTrait`](@ref)
"""
abstract type SubscribableTrait end

"""
Simple subscribable trait behavior. Simple subscribable can be used in subscribe! function and just executes `on_subscribe!` method for provided subscribable.
`SimpleSubscribableTrait` is a subtype of `SubscribableTrait`.

See also: [`SubscribableTrait`](@ref), [`Subscribable`](@ref), [`subscribe!`](@ref), [`on_subscribe!`](@ref)
"""
struct SimpleSubscribableTrait{T} <: SubscribableTrait end

"""
Scheduled subscribable trait behavior. Scheduled subscribable can be used in subscribe! function and executes `on_subscribe!` method for provided subscribable with custom scheduling.
`ScheduledSubscribableTrait` is a subtype of `SubscribableTrait`.

See also: [`SubscribableTrait`](@ref), [`ScheduledSubscribable`](@ref), [`subscribe!`](@ref), [`on_subscribe!`](@ref)
"""
struct ScheduledSubscribableTrait{T} <: SubscribableTrait end

"""
Default subscribable trait behavior for all types. Invalid subscribable cannot be used in subscribe! function, doing so will throw an error.
`InvalidSubscribableTrait` is a subtype of `SubscribableTrait`.

See also: [`SubscribableTrait`](@ref), [`subscribe!`](@ref)
"""
struct InvalidSubscribableTrait <: SubscribableTrait end

"""
Supertype type for `Subscribable` and `ScheduledSubscribable` types.

See also: [`Subscribable`](@ref), [`ScheduledSubscribable`](@ref)
"""
abstract type AbstractSubscribable{T} end

"""
Super type for any simple subscribable object. Automatically specifies a `SimpleSubscribableTrait` trait behavior.
Objects with specified `SimpleSubscribableTrait` subscribable trait must implement: `on_subscribe!(subscribable, actor)` method.
`Subscribable` is a subtype of `AbstractSubscribable` type.

# Examples
```jldoctest
using Rocket

struct MySubscribable <: Subscribable{String} end

Rocket.as_subscribable(MySubscribable)

# output

SimpleSubscribableTrait{String}()
```

See also: [`SubscribableTrait`](@ref), [`SimpleSubscribableTrait`](@ref)
"""
abstract type Subscribable{T} <: AbstractSubscribable{T} end

"""
Super type for any scheduled subscribable object. Automatically specifies a `ScheduledSubscribableTrait` trait behavior.
Objects with specified `ScheduledSubscribableTrait` subscribable trait must implement: `on_subscribe!(subscribable, actor, scheduler)` method.
`ScheduledSubscribable` is a subtype of `AbstractSubscribable` type.

# Examples
```jldoctest
using Rocket

struct MyScheduledSubscribable <: ScheduledSubscribable{String} end

Rocket.as_subscribable(MyScheduledSubscribable)

# output

ScheduledSubscribableTrait{String}()
```

See also: [`SubscribableTrait`](@ref), [`ScheduledSubscribableTrait`](@ref)
"""
abstract type ScheduledSubscribable{T} <: AbstractSubscribable{T} end

"""
    as_subscribable(any)

This function checks subscribable trait behavior specification. Can be used explicitly to specify subscribable trait behavior for any object.

See also: [`SubscribableTrait`](@ref)
"""
as_subscribable(::Type)                                        = InvalidSubscribableTrait()
as_subscribable(::Type{ <: Subscribable{T} })          where T = SimpleSubscribableTrait{T}()
as_subscribable(::Type{ <: ScheduledSubscribable{T} }) where T = ScheduledSubscribableTrait{T}()
as_subscribable(::O)                                   where O = as_subscribable(O)

subscribable_extract_type(::Type{S}) where S = subscribable_extract_type(as_subscribable(S), S)
subscribable_extract_type(::S)       where S = subscribable_extract_type(S)

subscribable_extract_type(::SimpleSubscribableTrait{T},    _) where T = T
subscribable_extract_type(::ScheduledSubscribableTrait{T}, _) where T = T
subscribable_extract_type(::InvalidSubscribableTrait, source)         = throw(InvalidSubscribableTraitUsageError(source))

Base.eltype(::Subscribable{T})          where T = T
Base.eltype(::ScheduledSubscribable{T}) where T = T

Base.eltype(::Type{ <: Subscribable{T} })          where T = T
Base.eltype(::Type{ <: ScheduledSubscribable{T} }) where T = T

union_type(::AbstractVector{ <: AbstractSubscribable{T} }) where T = T

"""
    subscribe!(subscribable::T, actor::S)   where { T, S }
    subscribe!(subscribable::T, factory::F) where { T, F <: AbstractActorFactory }
    subscribe!(subscriptions::Tuple)
    subscribe!(subscriptions::AbstractVector)

`subscribe!` function is used to attach an actor to subscribable.
It also checks types of subscribable and actors to be a valid Subscribable and Actor objects respectively.
Passing not valid subscribable or/and actor object will throw an error. If the input argument to the `subscribe!` function is either a tuple or a vector,
it will first check that all of the arguments are valid source and actor objects
and if its true will subscribe for each of them individually.

# Arguments
- `subscribable`: valid subscribable object
- `actor`: valid actor object

# Examples

```jldoctest
using Rocket

source = from((1, 2, 3))
subscribe!(source, logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed
```

```jldoctest
using Rocket

source = from((1, 2, 3))
subscribe!(source, 1)
;

# output

ERROR: Type Int64 is not a valid actor type.
[...]
```

```jldoctest
using Rocket

source = from((1, 2, 3))
subscribe!(1, logger())
;

# output

ERROR: Type Int64 is not a valid subscribable type.
[...]
```

See also: [`on_subscribe!`](@ref), [`as_subscribable`](@ref)
"""
@inline subscribe!(subscribable::T, actor::S) where { T, S } = check_on_subscribe!(as_subscribable(T), as_actor(S), subscribable, actor)

# Specialised methods for built-in default subscribable and actor types
@inline subscribe!(subscribable::Subscribable{T1},          actor::Actor{T2}) where { T2, T1 <: T2 } = on_subscribe!(subscribable, actor)
@inline subscribe!(subscribable::ScheduledSubscribable{T1}, actor::Actor{T2}) where { T2, T1 <: T2 } = scheduled_subscription!(subscribable, actor, makeinstance(T1, getscheduler(subscribable)))

@inline subscribe!(subscribable::Subscribable{T},          actor::Actor{T}) where { T } = on_subscribe!(subscribable, actor)
@inline subscribe!(subscribable::ScheduledSubscribable{T}, actor::Actor{T}) where { T } = scheduled_subscription!(subscribable, actor, makeinstance(T, getscheduler(subscribable)))

function subscribe!(args::Union{Tuple, AbstractVector})
    # First check that all arguments are correct
    foreach(enumerate(args)) do (index, arg)
        if length(arg) === 2
            source = first(arg)
            actor  = last(arg)
            if as_subscribable(typeof(source)) === InvalidSubscribableTrait()
                error("multiple_subscribe! usage: Invalid source '$(source)' found on index $(index)")
            elseif !(typeof(actor) <: AbstractActorFactory) && !(typeof(actor) <: Function) && as_actor(typeof(actor)) === InvalidActorTrait()
                error("multiple_subscribe! usage: Invalid actor '$(actor)' found on index $(index)")
            end
        else
            error("Invalid multiple_subscribe! usage: Individual argument should be a size-2 Tuple, found '$(arg)' on index $(index)")
        end
    end
    return map(args) do arg
        try
            return subscribe!(first(arg), last(arg))
        catch error
            @error "Error occured during multiple subscription. Return `voidTeardown` for '$arg' argument."
            @error error
            return voidTeardown
        end
    end
end

# We don't use an abstract types here and dispatch on all possible combinations of types because of the issue #37045 JuliaLang/julia
# https://github.com/JuliaLang/julia/issues/37045
@inline check_on_subscribe!(::InvalidSubscribableTrait,       _,                          subscribable, actor)                  = throw(InvalidSubscribableTraitUsageError(subscribable))
@inline check_on_subscribe!(::SimpleSubscribableTrait,        ::InvalidActorTrait,        subscribable, actor)                  = throw(InvalidActorTraitUsageError(actor))
@inline check_on_subscribe!(::ScheduledSubscribableTrait,     ::InvalidActorTrait,        subscribable, actor)                  = throw(InvalidActorTraitUsageError(actor))
@inline check_on_subscribe!(::SimpleSubscribableTrait{T1},    ::BaseActorTrait{T2},       subscribable, actor) where { T1, T2 } = throw(InconsistentActorWithSubscribableDataTypesError{T1, T2}(subscribable, actor))
@inline check_on_subscribe!(::SimpleSubscribableTrait{T1},    ::NextActorTrait{T2},       subscribable, actor) where { T1, T2 } = throw(InconsistentActorWithSubscribableDataTypesError{T1, T2}(subscribable, actor))
@inline check_on_subscribe!(::SimpleSubscribableTrait{T1},    ::ErrorActorTrait{T2},      subscribable, actor) where { T1, T2 } = throw(InconsistentActorWithSubscribableDataTypesError{T1, T2}(subscribable, actor))
@inline check_on_subscribe!(::SimpleSubscribableTrait{T1},    ::CompletionActorTrait{T2}, subscribable, actor) where { T1, T2 } = throw(InconsistentActorWithSubscribableDataTypesError{T1, T2}(subscribable, actor))
@inline check_on_subscribe!(::ScheduledSubscribableTrait{T1}, ::BaseActorTrait{T2},       subscribable, actor) where { T1, T2 } = throw(InconsistentActorWithSubscribableDataTypesError{T1, T2}(subscribable, actor))
@inline check_on_subscribe!(::ScheduledSubscribableTrait{T1}, ::NextActorTrait{T2},       subscribable, actor) where { T1, T2 } = throw(InconsistentActorWithSubscribableDataTypesError{T1, T2}(subscribable, actor))
@inline check_on_subscribe!(::ScheduledSubscribableTrait{T1}, ::ErrorActorTrait{T2},      subscribable, actor) where { T1, T2 } = throw(InconsistentActorWithSubscribableDataTypesError{T1, T2}(subscribable, actor))
@inline check_on_subscribe!(::ScheduledSubscribableTrait{T1}, ::CompletionActorTrait{T2}, subscribable, actor) where { T1, T2 } = throw(InconsistentActorWithSubscribableDataTypesError{T1, T2}(subscribable, actor))

@inline check_on_subscribe!(::SimpleSubscribableTrait{T1}, ::BaseActorTrait{T2},       subscribable, actor) where { T2, T1 <: T2 } = on_subscribe!(subscribable, actor)
@inline check_on_subscribe!(::SimpleSubscribableTrait{T1}, ::NextActorTrait{T2},       subscribable, actor) where { T2, T1 <: T2 } = on_subscribe!(subscribable, actor)
@inline check_on_subscribe!(::SimpleSubscribableTrait{T1}, ::ErrorActorTrait{T2},      subscribable, actor) where { T2, T1 <: T2 } = on_subscribe!(subscribable, actor)
@inline check_on_subscribe!(::SimpleSubscribableTrait{T1}, ::CompletionActorTrait{T2}, subscribable, actor) where { T2, T1 <: T2 } = on_subscribe!(subscribable, actor)

@inline check_on_subscribe!(::ScheduledSubscribableTrait{T1}, ::BaseActorTrait{T2},       subscribable, actor) where { T2, T1 <: T2 } = scheduled_subscription!(subscribable, actor, makeinstance(T1, getscheduler(subscribable)))
@inline check_on_subscribe!(::ScheduledSubscribableTrait{T1}, ::NextActorTrait{T2},       subscribable, actor) where { T2, T1 <: T2 } = scheduled_subscription!(subscribable, actor, makeinstance(T1, getscheduler(subscribable)))
@inline check_on_subscribe!(::ScheduledSubscribableTrait{T1}, ::ErrorActorTrait{T2},      subscribable, actor) where { T2, T1 <: T2 } = scheduled_subscription!(subscribable, actor, makeinstance(T1, getscheduler(subscribable)))
@inline check_on_subscribe!(::ScheduledSubscribableTrait{T1}, ::CompletionActorTrait{T2}, subscribable, actor) where { T2, T1 <: T2 } = scheduled_subscription!(subscribable, actor, makeinstance(T1, getscheduler(subscribable)))

@inline function subscribe!(subscribable::T, factory::F) where { T, F <: AbstractActorFactory }
    return check_on_subscribe_with_factory!(as_subscribable(T), subscribable, factory)
end

@inline check_on_subscribe_with_factory!(::InvalidSubscribableTrait,      subscribable, factory)         = throw(InvalidSubscribableTraitUsageError(subscribable))
@inline check_on_subscribe_with_factory!(::SimpleSubscribableTrait{L},    subscribable, factory) where L = subscribe!(subscribable, create_actor(L, factory))
@inline check_on_subscribe_with_factory!(::ScheduledSubscribableTrait{L}, subscribable, factory) where L = subscribe!(subscribable, create_actor(L, factory))

"""
    on_subscribe!(subscribable, actor)
    on_subscribe!(subscribable, actor, scheduler)

Every valid subscribable object have to define its own method for `on_subscribe!` function which specifies subscription logic
and has return a valid `Teardown` object.

Objects with specified `SimpleSubscribableTrait` subscribable trait must implement: `on_subscribe!(subscribable, actor)` method.
Objects with specified `ScheduledSubscribableTrait` subscribable trait must implement: `on_subscribe!(subscribable, actor, scheduler)` method.

# Arguments
- `subscribable`: Subscribable object
- `actor`: Actor object
- `scheduler`: Scheduler object (only for scheduled subscribables)

# Examples

```jldoctest
using Rocket

struct MySubscribable <: Subscribable{Int} end

function Rocket.on_subscribe!(subscribable::MySubscribable, actor)
    next!(actor, 0)
    complete!(actor)
    return voidTeardown
end

subscribe!(MySubscribable(), logger())
;

# output

[LogActor] Data: 0
[LogActor] Completed
```

```jldoctest
using Rocket

struct MyScheduledSubscribable <: ScheduledSubscribable{Int} end

Rocket.getscheduler(::MyScheduledSubscribable) = AsapScheduler()

function Rocket.on_subscribe!(subscribable::MyScheduledSubscribable, actor, scheduler)
    next!(actor, 0, scheduler)
    complete!(actor, scheduler)
    return voidTeardown
end

subscribe!(MyScheduledSubscribable(), logger())
;

# output

[LogActor] Data: 0
[LogActor] Completed
```

See also: [`Subscribable`](@ref), [`ScheduledSubscribable`](@ref), [`SimpleSubscribableTrait`](@ref), [`ScheduledSubscribableTrait`](@ref), [`Teardown`](@ref), [`logger`](@ref)
"""
on_subscribe!(subscribable, actor)            = throw(MissingOnSubscribeImplementationError(subscribable))
on_subscribe!(subscribable, actor, scheduler) = throw(MissingOnScheduledSubscribeImplementationError(subscribable))

# -------------------------------- #
# Errors                           #
# -------------------------------- #

"""
This error will be thrown if `subscribe!` function is called with invalid subscribable object

See also: [`subscribe!`](@ref)
"""
struct InvalidSubscribableTraitUsageError
    subscribable
end

function Base.showerror(io::IO, err::InvalidSubscribableTraitUsageError)
    print(io, "Type $(typeof(err.subscribable)) is not a valid subscribable type. \nConsider extending your subscribable with Subscribable{T} abstract type or implement as_subscribable(::Type{<:$(typeof(err.subscribable))}).")
end

"""
This error will be thrown if `subscribe!` function is called with inconsistent subscribable and actor objects

See also: [`subscribe!`](@ref)
"""
struct InconsistentActorWithSubscribableDataTypesError{T1, T2}
    subscribable
    actor
end

function Base.showerror(io::IO, err::InconsistentActorWithSubscribableDataTypesError{T1, T2}) where T1 where T2
    print(io, "Actor of type $(typeof(err.actor)) expects data to be of type $(T2), while subscribable of type $(typeof(err.subscribable)) produces data of type $(T1).")
end

"""
This error will be thrown if Julia cannot find specific method of 'on_subscribe!()' function for given subscribable and actor

See also: [`on_subscribe!`](@ref)
"""
struct MissingOnSubscribeImplementationError
    subscribable
end

function Base.showerror(io::IO, err::MissingOnSubscribeImplementationError)
    print(io, "You probably forgot to implement on_subscribe!(subscribable::$(typeof(err.subscribable)), actor).")
end

"""
This error will be thrown if Julia cannot find specific method of 'on_subscribe!()' function for given subscribable, actor and scheduler

See also: [`on_subscribe!`](@ref)
"""
struct MissingOnScheduledSubscribeImplementationError
    subscribable
end

function Base.showerror(io::IO, err::MissingOnScheduledSubscribeImplementationError)
    print(io, "You probably forgot to implement on_subscribe!(subscribable::$(typeof(err.subscribable)), actor, scheduler).")
end


text:
using LightGraphs
using GraphBenchmarks
import GraphBenchmarks: construct, preparebench
using BenchmarkTools

export LGGraph, construct

immutable LGGraph{T<:SimpleGraph} <: GraphType
    g::T
end

function construct{T <: SimpleGraph}(t::Type{LGGraph{T}}, edges::Array{Int64, 2})
    g = T(size(edges, 2))
    for i=1:size(edges,2)
        add_edge!(g, edges[1, i], edges[2, i])
    end
    LGGraph(g)
end

function preparebench(benchmark::GraphBenchmarkSpec, alg::BFS, lg::LGGraph)
    n = nv(lg.g)
    sources = picksources(benchmark, n)
    visitor = LightGraphs.TreeBFSVisitorVector(zeros(Int,n))
    bfsbench = @benchmarkable begin
        for src in $sources
            LightGraphs.bfs_tree!($visitor, $(lg.g), src)
        end
    end seconds=6000 samples=10
    bfsbench
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5061fe80-bbdf-443c-ab81-ac1a6e2d5160
using BSON, Flux, ProgressMeter, Setfield, Random, ROC, NPZ

# ╔═╡ f375b5a5-cb41-40c1-be03-e413e744c6be
using GraphFlux

# ╔═╡ 93e83138-d84a-11eb-1c94-c54c588f853d
# begin
# 	using Pkg

# 	Pkg.add(["BSON", "Flux", "ProgressMeter", "Setfield", "Random", "NPZ"])
# 	Pkg.add(url="https://github.com/diegozea/ROC.jl")
# end

# ╔═╡ e20d7b3d-6559-477d-a258-7fcfe690837b
md"""If you need to download the data the first time, uncomment the hidden cell below"""

# ╔═╡ f2e46818-2bce-4a70-9908-933b638e42d0
# begin
# 	run(`curl http://snap.stanford.edu/ogb/data/graphproppred/csv_mol_download/hiv.zip --create-dirs -o "../download/hiv.zip"`)
# 	run(`unzip ../download/hiv.zip -d ../download`)

# 	Pkg.add(["CSV", "CodecZlib", "Mmap", "DataFrames"])
# 	using CSV, CodecZlib, Mmap, DataFrames

# 	function csv_to_arr(fname)
# 		memmap = Mmap.mmap(fname)
# 		csv = CSV.File(transcode(GzipDecompressor, memmap), header=false)
# 		csv |> DataFrame |> Array
# 	end

# 	prefix = "../download/hiv/raw/"

# 	edgelist = csv_to_arr(prefix * "edge.csv.gz") .+ 1

# 	nodes = csv_to_arr(prefix * "node-feat.csv.gz")' .+ 1
# 	edges = csv_to_arr(prefix * "edge-feat.csv.gz")' .+ 1

# 	trainidx = csv_to_arr("../download/hiv/split/scaffold/train.csv.gz")[:, 1] .+ 1
# 	validx = csv_to_arr("../download/hiv/split/scaffold/valid.csv.gz")[:, 1] .+ 1
# 	testidx = csv_to_arr("../download/hiv/split/scaffold/test.csv.gz")[:, 1] .+ 1

# 	targets = csv_to_arr(prefix * "graph-label.csv.gz")[:, 1]

# 	numnodes = csv_to_arr(prefix * "num-node-list.csv.gz")[:, 1]
# 	numedges = csv_to_arr(prefix * "num-edge-list.csv.gz")[:, 1]

# 	graphdata = Dict("numnodes" => numnodes, "numedges" => numedges, "edgelist" => edgelist)

# 	features = Dict("nodes" => nodes, "edges" => edges, "targets" => targets)
# 	idxs = Dict("train" => trainidx, "valid" => validx, "test" => testidx)

# 	BSON.@save "../download/hiv.bson" graphdata features idxs
# end

# ╔═╡ 5adf7e4f-0eee-482c-a6d5-fb85b0de7bc3
BSON.@load "../download/hiv.bson" features graphdata idxs

# ╔═╡ 238a3c8d-70cd-4077-8804-68a5449da136
md"""We define the model here, leveraging heavily from `features/molecules/ogbhiv.jl`"""

# ╔═╡ 07cf68c9-6e8d-4316-92aa-df879b690039
begin
	struct OGBGCNBlock
	    gcn
	    batchnorm
	    dropout
	    last
	end

	function OGBGCNBlock(embdim::Integer; droprate::AbstractFloat=0.5, last::Bool=false)
	    OGBGCNBlock(GCNₑ(embdim, embdim),
	                BatchNorm(embdim),
	                Dropout(droprate),
	                last)
	end

	function (block::OGBGCNBlock)(g)
	    g = block.gcn(g)
	    nodeh = block.batchnorm(GraphFlux.nodes(g))
	    if block.last nodeh = relu.(nodeh) end
	    nodeh = block.dropout(nodeh)
	    @set g.nodes = nodeh
	end

	Flux.@functor OGBGCNBlock
end

# ╔═╡ dda37dd5-b122-4666-acf2-a7750863fd2e
# model constructor
model() = Chain(GraphFlux.OGBAtomEncoder(300),
	            OGBGCNBlock(300),
	            OGBGCNBlock(300),
	            OGBGCNBlock(300),
	            OGBGCNBlock(300),
	            OGBGCNBlock(300; last=true),
	            graphmeanpool,
	            Dense(300, 1),
	            x -> x[1, :])

# ╔═╡ 5dd07e15-f9a5-41f4-ae79-b5eee2ba0087
md"""Total number of parameters match that used in the original paper :)"""

# ╔═╡ 3a9e78d2-5514-48b9-a1cf-3370fe92ebae
# number of parameters in model, total
[length(x) for x in Flux.params(model())] |> sum

# ╔═╡ 72bace50-0338-4f95-8c43-f74bf45d236e
md"""We define a training run here. There's probably a cleaner way of doing this, but this was fast."""

# ╔═╡ 947b7429-efb5-4cc5-a5d8-e72ba5c7fb97
begin
	NUM_TRIALS = 5
	EPOCHS = 100
	TRAIN_BATCHSIZE = 1024
end

# ╔═╡ a26bcee2-5c06-452a-bb79-9382890b5404
md"""We define temporary evaluation functions (e.g. `getroc`) in the hidden cell here, final evaluations are done inside Python with OGB Evaluators"""

# ╔═╡ ca67a788-5c58-42bf-8352-8ad7f3984b12
begin
	function getpreds(m, subset="valid")
		vcat([m(GraphFlux.getmolbatch(features, graphdata, idsub) |> gpu) for idsub in Iterators.partition(idxs[subset], 2048)]...)
	end

	function getroc(m, subset="valid")
		preds = getpreds(m, subset)
		AUC(roc(sigmoid.(preds) |> cpu, features["targets"][idxs[subset]]))
	end
end

# ╔═╡ 8cb4646b-bce3-4d21-a50c-d34cdcb6b64a
function training_run(model_num)
	# create a new model
	m = model() |> gpu
	ps = Flux.params(m)

	# Metrics need to be computed on model with best val roc as in here:
	# https://github.com/snap-stanford/ogb/blob/master/examples/graphproppred/mol/main_pyg.py
	bestvalroc = 0
	bestmodel = nothing

	# new loss fn and optimizers
	loss(x, y) = Flux.Losses.logitbinarycrossentropy(m(x), y)
	opt = Flux.Optimise.ADAM()

	allbatchidxs = Iterators.partition(Random.shuffle(idxs["train"]), TRAIN_BATCHSIZE)

	for epoch in 1:EPOCHS
		@showprogress for batchidxs in allbatchidxs
			X = GraphFlux.getmolbatch(features, graphdata, batchidxs) |> gpu
			y = features["targets"][batchidxs] |> gpu

			data = [(X, y)]
			Flux.train!(loss, ps, data, opt)
		end

		# visible in terminal
		rocauc = getroc(m, "valid")
		@show (epoch, rocauc)

		if rocauc > bestvalroc
			@show bestvalroc = rocauc
			bestmodel = deepcopy(m)
		end
	end

	finishedmodel = bestmodel |> cpu
	BSON.@save "../download/modelhiv$(model_num).bson" finishedmodel
	testpreds = sigmoid.(getpreds(bestmodel, "test")) |> cpu
	valpreds = sigmoid.(getpreds(bestmodel, "valid")) |> cpu

	return valpreds, testpreds
end

# ╔═╡ c4871692-9a1f-40c6-b6f2-5d069a753918
begin
	allvalpreds = []
	alltestpreds = []

	for i in 1:NUM_TRIALS
		@show "training run $(i)"
		valpreds, testpreds = training_run(i)
		push!(alltestpreds, testpreds)
		push!(allvalpreds, valpreds)
	end
end

# ╔═╡ 3b9abe03-f709-443e-b463-1fe839e5a65b
npzwrite("../download/hivpreds.npz",
	Dict("testpreds" => hcat(alltestpreds...),
		"testtargets" => features["targets"][idxs["test"]],
		"valpreds" => hcat(allvalpreds...),
		"valtargets" => features["targets"][idxs["valid"]]))

# ╔═╡ Cell order:
# ╠═93e83138-d84a-11eb-1c94-c54c588f853d
# ╠═5061fe80-bbdf-443c-ab81-ac1a6e2d5160
# ╟─e20d7b3d-6559-477d-a258-7fcfe690837b
# ╟─f2e46818-2bce-4a70-9908-933b638e42d0
# ╠═f375b5a5-cb41-40c1-be03-e413e744c6be
# ╠═5adf7e4f-0eee-482c-a6d5-fb85b0de7bc3
# ╟─238a3c8d-70cd-4077-8804-68a5449da136
# ╠═07cf68c9-6e8d-4316-92aa-df879b690039
# ╠═dda37dd5-b122-4666-acf2-a7750863fd2e
# ╟─5dd07e15-f9a5-41f4-ae79-b5eee2ba0087
# ╠═3a9e78d2-5514-48b9-a1cf-3370fe92ebae
# ╟─72bace50-0338-4f95-8c43-f74bf45d236e
# ╠═947b7429-efb5-4cc5-a5d8-e72ba5c7fb97
# ╠═8cb4646b-bce3-4d21-a50c-d34cdcb6b64a
# ╟─a26bcee2-5c06-452a-bb79-9382890b5404
# ╟─ca67a788-5c58-42bf-8352-8ad7f3984b12
# ╠═c4871692-9a1f-40c6-b6f2-5d069a753918
# ╠═3b9abe03-f709-443e-b463-1fe839e5a65b


text:
for PT in (:Particles, :StaticParticles)
    @eval begin
        Base.promote_rule(::Type{Complex{S}}, ::Type{$PT{T,N}}) where {S<:Real,T<:Real,N} = Complex{$PT{promote_type(S,T),N}}

    end

    for ff in (^,)
        f = nameof(ff)
        @eval Base.$f(z::Complex{$PT{T,N}}, x::Real) where {T,N} = ℂ2ℂ_function($f, z, x)
        @eval Base.$f(z::Real, x::Complex{$PT{T,N}}) where {T,N} = ℂ2ℂ_function($f, z, x)
        @eval Base.$f(z::Complex{$PT{T,N}}, x::Complex{$PT{T,N}}) where {T,N} = ℂ2ℂ_function($f, z, x)

        @eval Base.$f(z::Complex{$PT{T,N}}, x::Int) where {T,N} = ℂ2ℂ_function($f, z, x)
        @eval Base.$f(z::Int, x::Complex{$PT{T,N}}) where {T,N} = ℂ2ℂ_function($f, z, x)
    end

end

@inline maybe_complex_particles(x,i) = x
@inline maybe_complex_particles(p::AbstractParticles,i) = complex(real(p.particles[i]), imag(p.particles[i]))

"""
    ℂ2ℂ_function(f::Function, z::Complex{<:AbstractParticles})

Helper function for uncertainty propagation through complex-valued functions of complex arguments.
applies `f : ℂ → ℂ ` to `z::Complex{<:AbstractParticles}`.
"""
function ℂ2ℂ_function(f::F, z::Complex{T}) where {F<:Union{Function,DataType},T<:AbstractParticles}
    rz,iz = z.re,z.im
    s = map(1:length(rz.particles)) do i
        @inbounds f(complex(rz.particles[i], iz.particles[i]))
    end
    complex(T(real.(s)), T(imag.(s)))
end

function ℂ2ℂ_function(f::F, z::Union{Complex{T},T}, a::R) where {F<:Union{Function,DataType},T<:AbstractParticles,R<:Real}
    rz,iz = z.re,z.im
    s = map(1:length(rz.particles)) do i
        @inbounds f(complex(rz.particles[i], iz.particles[i]),  a)
    end
    complex(T(real.(s)), T(imag.(s)))
end

function ℂ2ℂ_function(f::F, z::R, a::Union{Complex{S},S}) where {F<:Union{Function,DataType},S<:AbstractParticles,R<:Real}
    rz,iz = a.re,a.im
    s = map(1:length(rz.particles)) do i
        @inbounds f(z, complex(rz.particles[i], iz.particles[i]))
    end
    complex(S(real.(s)), S(imag.(s)))
end

function ℂ2ℂ_function(f::F, z::Union{Complex{T},T}, a::Union{Complex{S},S}) where {F<:Union{Function,DataType},T<:AbstractParticles,S<:AbstractParticles}
    rz,iz = z.re,z.im
    ra,ia = a.re,a.im

    s = map(1:length(rz.particles)) do i
        @inbounds f(complex(rz.particles[i], iz.particles[i]),  complex(ra.particles[i], ia.particles[i]))
    end
    complex(T(real.(s)), T(imag.(s)))
end

function ℂ2ℂ_function!(f::F, s, z::Complex{T}) where {F,T<:AbstractParticles}
    rz,iz = z.re,z.im
    map!(s, 1:length(rz.particles)) do i
        @inbounds f(complex(rz.particles[i], iz.particles[i]))
    end
    complex(T(real.(s)), T(imag.(s)))
end

for ff in (sqrt, exp, exp10, log, log10, sin, cos, tan)
    f = nameof(ff)
    @eval Base.$f(z::Complex{<: AbstractParticles}) = ℂ2ℂ_function($f, z)
    @eval $(Symbol(f,:!))(s, z::Complex{<: AbstractParticles}) = ℂ2ℂ_function!($f, s, z)
end

Base.isinf(p::Complex{<: AbstractParticles}) = isinf(real(p)) || isinf(imag(p))
Base.isfinite(p::Complex{<: AbstractParticles}) = isfinite(real(p)) && isfinite(imag(p))

function Base.:(/)(a::Complex{T}, b::Complex{T}) where T<:AbstractParticles
    are = real(a); aim = imag(a); bre = real(b); bim = imag(b)
    if mean(abs(bre)) <= mean(abs(bim))
        if isinf(bre) && isinf(bim)
            r = sign(bre)/sign(bim)
        else
            r = bre / bim
        end
        den = bim + r*bre
        Complex((are*r + aim)/den, (aim*r - are)/den)
    else
        if isinf(bre) && isinf(bim)
            r = sign(bim)/sign(bre)
        else
            r = bim / bre
        end
        den = bre + r*bim
        Complex((are + aim*r)/den, (aim - are*r)/den)
    end
end

function Base.:(/)(are::T, b::Complex{T}) where T<:AbstractParticles
    aim = 0; bre = real(b); bim = imag(b)
    if mean(abs(bre)) <= mean(abs(bim))
        if isinf(bre) && isinf(bim)
            r = sign(bre)/sign(bim)
        else
            r = bre / bim
        end
        den = bim + r*bre
        Complex((are*r)/den, (-are)/den)
    else
        if isinf(bre) && isinf(bim)
            r = sign(bim)/sign(bre)
        else
            r = bim / bre
        end
        den = bre + r*bim
        Complex((are)/den, (-are*r)/den)
    end
end


text:
include("help.jl")
function mark_angles(r)
    num_steps=[]
    while isborder(r,Sud)==false || isborder(r,West) == false # Робот - не в юго-западном углу
        push!(num_steps, moves!(r, West))
        push!(num_steps, moves!(r, Sud))
    end
    #Робот - в юго-западном углу и в num_steps - закодирован пройденный путь
    for side in (Nord,Ost,Sud,West)
        movements!(r,side)
        putmarker!(r)
    end
    # Маркеры поставлены и Робот - в юго-западном углу

    for (i,n) in enumerate(num_steps)
        side = isodd(i) ? Ost : Nord
        movements!(r,side,n)
    end
    # Робот - в исходном положении
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
module Cards

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element

export card, card_section, card_actions

function __init__()
  Genie.Renderer.Html.register_normal_element("q__card", context = Genie.Renderer.Html)
  Genie.Renderer.Html.register_normal_element("q__card__section", context = Genie.Renderer.Html)
  Genie.Renderer.Html.register_normal_element("q__card__actions", context = Genie.Renderer.Html)
end

function card(args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  wrap() do
    Genie.Renderer.Html.q__card(args...; kwargs...)
  end
end

function card_section(args...; kwargs...)
  Genie.Renderer.Html.q__card__section(args...; kwargs...)
end

function card_actions(args...; kwargs...)
  Genie.Renderer.Html.q__card__actions(args...; kwargs...)
end

end


text:
function create_pgmfile(img, outf::String)
    s = open(outf, "w")
    write(s, "P5\n")
    n, m = size(img)
    write(s, "$m $n 255\n")
    for i=1:n, j=1:m
        p = img[i,j]
        write(s, uint8(p))
    end
    close(s)
end

max_stars_repo_path: Module 3/Chapter02/pgmfile.jl

max_issues_repo_path: Module 3/Chapter02/pgmfile.jl

max_forks_repo_path: Module 3/Chapter02/pgmfile.jl


text:
function tgear(thtl)  # Power command v. thtl. relationship
    if thtl <= 0.77
        cpow = 64.94 * thtl
    else
        cpow = 217.38 * thtl - 117.38
    end
    return cpow
end

function pdot(p3, p1) # PDOT= rate of change of power
    if p1 >= 50.0  # P3= actual power, P1= power command
        if p3 >= 50.0
            t = 5.0
            p2 = p1
        else
            p2 = 60.0
            t = rtau(p2 - p3)
        end
    else
        if p3 >= 50.0
            t = 5.0
            p2 = 40.0
        else
            p2 = p1
            t = rtau(p2 - p3)
        end
    end
    pd = t * (p2 - p3)
    return pd
end

function rtau(dp)
    if dp <= 25.0
        rt = 1.0  # Reciprocal time constant
    elseif dp >= 50.0
        rt = 0.1
    else
        rt = 1.9 - 0.036 * dp
    end
    return rt
end

idle_data =
    [
        1060.0 670.0 880.0 1140.0 1500.0 1860.0
        635.0 425.0 690.0 1010.0 1330.0 1700.0
        60.0 25.0 345.0 755.0 1130.0 1525.0
        -1020.0 -710.0 -300.0 350.0 910.0 1360.0
        -2700.0 -1900.0 -1300.0 -247.0 600.0 1100.0
        -3600.0 -1400.0 -595.0 -342.0 -200.0 700.0
    ]'

mil_data =
    [
        12680.0 9150.0 6200.0 3950.0 2450.0 1400.0
        12680.0 9150.0 6313.0 4040.0 2470.0 1400.0
        12610.0 9312.0 6610.0 4290.0 2600.0 1560.0
        12640.0 9839.0 7090.0 4660.0 2840.0 1660.0
        12390.0 10176.0 7750.0 5320.0 3250.0 1930.0
        11680.0 9848.0 8050.0 6100.0 3800.0 2310.0
    ]'

max_data =
    [
        20000.0 15000.0 10800.0 7000.0 4000.0 2500.0
        21420.0 15700.0 11225.0 7323.0 4435.0 2600.0
        22700.0 16860.0 12250.0 8154.0 5000.0 2835.0
        24240.0 18910.0 13760.0 9285.0 5700.0 3215.0
        26070.0 21075.0 15975.0 11115.0 6860.0 3950.0
        28886.0 23319.0 18300.0 13484.0 8642.0 5057.0
    ]'

function thrust(pow, alt, rmach)  # Engine thrust model

    alt = alt * M2FT

    if alt <= 0.
        alt = 0.0
    end

    h = 0.0001 * alt
    i = floor(Int, h)
    if i >= 5
        i = 4
    end

    dh = h - float(i)

    rm = 5.0 * rmach
    m = floor(Int, rm)
    if m >= 5
        m = 4
    end

    dm = rm - float(m)
    cdh = 1.0 - float(dh)

    i = i + 1
    m = m + 1

    s = mil_data[i, m] * cdh + mil_data[i+1, m] * dh
    t = mil_data[i, m+1] * cdh + mil_data[i+1, m+1] * dh
    tmil = s + (t - s) * dm

    if pow < 50.0
        s = idle_data[i, m] * cdh + idle_data[i+1, m] * dh
        t = idle_data[i, m+1] * cdh + idle_data[i+1, m+1] * dh
        tidl = s + (t - s) * dm
        thrst = tidl + (tmil - tidl) * pow * 0.02
    else
        s = max_data[i, m] * cdh + max_data[i+1, m] * dh
        t = max_data[i, m+1] * cdh + max_data[i+1, m+1] * dh
        tmax = s + (t - s) * dm
        thrst = tmil + (tmax - tmil) * (pow - 50.0) * 0.02
    end

    # Transform from lbf to N (using GD defined in Stevens)
    thrst = thrst * LB2KG * GD * FT2M
    return thrst
end

function calculate_prop_forces_moments(ac::F16, dss::DSState, mach, controls)

    # Assign state & control variables
    alt = get_height(dss)
    pow = get_engine_power(dss)

    thtl = controls[1]

    thrust_ = thrust(pow, alt, mach)

    return [thrust_, 0.0, 0.0, 0.0, 0.0, 0.0]
end

calculate_prop_gyro_effects(ac::F16) = [HX, 0, 0]

calculate_pdot(ac::F16, thtl, pow) = pdot(pow, tgear(thtl))

tgear(ac::F16, thtl) = tgear(thtl)


text:
using Parameters, Compat

@with_kw mutable struct lsoda_common_t
    yh::Ptr{Ptr{Cdouble}}= C_NULL
    wm::Ptr{Ptr{Cdouble}}= C_NULL
    ewt::Ptr{Cdouble}= C_NULL
    savf::Ptr{Cdouble}= C_NULL
    acor::Ptr{Cdouble}= C_NULL
    ipvt::Ptr{Cint}= C_NULL
    memory::Ptr{Nothing}= C_NULL
    h::Cdouble = 0.
    hu::Cdouble = 0.
    rc::Cdouble = 0.
    tn::Cdouble = 0.
    tsw::Cdouble = 0.
    pdnorm::Cdouble = 0.
    crate::Cdouble = 0.
    el::NTuple{14,Cdouble} = ntuple(x->Cdouble(0), 14)
    elco::NTuple{13,NTuple{14,Cdouble}} = ntuple(x->ntuple(x->Cdouble(0), 14),13)
    tesco::NTuple{13,NTuple{4,Cdouble}} = ntuple(x->ntuple(x->Cdouble(0), 4),13)
    hold::Cdouble = 0.
    rmax::Cdouble = 0.
    pdest::Cdouble = 0.
    pdlast::Cdouble = 0.
    ialth::Cint = 0
    ipup::Cint = 0
    nslp::Cint = 0
    icount::Cint = 0
    irflag::Cint = 0
    imxer::Cint = 0
    illin::Cint = 0
    nhnil::Cint = 0
    nslast::Cint = 0
    jcur::Cint = 0
    meth::Cint = 0
    mused::Cint = 0
    nq::Cint = 0
    nst::Cint = 0
    ncf::Cint = 0
    nfe::Cint = 0
    nje::Cint = 0
    nqu::Cint = 0
    miter::Cint = 0
end

@with_kw mutable struct lsoda_opt_t
	  ixpr::Cint = 0
    mxstep::Cint = 0
    mxhnil::Cint = 0
    mxordn::Cint = 0
    mxords::Cint = 0
    tcrit::Cdouble = 0.
    h0::Cdouble = 0.
    hmax::Cdouble = 0.
    hmin::Cdouble = 0.
    hmxi::Cdouble = 0.
    itask::Cint = 0
    rtol::Ptr{Cdouble} = C_NULL
    atol::Ptr{Cdouble} = C_NULL
end

const _lsoda_f = Ptr{Nothing}

@with_kw mutable struct lsoda_context_t
    function_::_lsoda_f = C_NULL
    data::Ptr{Nothing} = C_NULL ##
    neq::Cint = 0
    state::Cint = 0
    error::Cstring = C_NULL
    common::Ptr{lsoda_common_t} = C_NULL
    opt::Ptr{lsoda_opt_t} = C_NULL
end

const lsoda_context_t_ptr = Ptr{lsoda_context_t}

mutable struct UserFunctionAndData
    func::Function
    data::Any
    neq::Cint
    UserFunctionAndData(func::Function, data::Any, neq::Cint) = new(func, data, neq)
end

# UserFunctionAndData(func::Function) = func
# UserFunctionAndData(func::Function, data::Nothing) = func
# UserFunctionAndData(func::Function, data::Nothing, neq::Cint) = func

function lsoda_prepare(ctx::lsoda_context_t,opt::lsoda_opt_t)
  return ccall((:lsoda_prepare,liblsoda),Cint,
    (Ptr{lsoda_context_t},Ptr{lsoda_opt_t}),
    Ref(ctx),Ref(opt))
end

function lsoda(ctx::lsoda_context_t,y::Vector,t::Vector{Float64},tout)
  return ccall((:lsoda,liblsoda),Cint,
    (Ptr{lsoda_context_t},Ptr{Cdouble},Ptr{Cdouble},Cdouble),
    Ref(ctx),y,t,tout[1])
end

function lsoda_reset(ctx::lsoda_context_t)
	ccall((:lsoda_reset,liblsoda),Nothing,(Ptr{lsoda_context_t},),Ref(ctx))
end

# written to wrap lsoda_free from C library but never used in practise as
# lsoda_context_t variables are handled on Julia's side
function lsoda_free(ctx::lsoda_context_t)
		ccall((:lsoda_free,liblsoda),Nothing,(Ref{lsoda_context_t},),Ref(ctx))
    nothing
end

function lsoda_free(ctx::Ref{lsoda_context_t})
		ccall((:lsoda_free,liblsoda),Nothing,(Ref{lsoda_context_t},),(ctx))
    nothing
end


text:
using Pkg
Pkg.add(url="https://github.com/LilithHafner/FBD")

using FBDCompare
using Test

ci = ("CI"=>"true") ∈ ENV

@testset "HyperPA" begin
    HyperPA.jl_benchmark()
    HyperPA.jl_profile()
    HyperPA.jl_with_io()
    if ci
        @test_broken false #=
https://github.com/LilithHafner/FBDCompare.jl/runs/4221035425?check_suite_focus=true
Traceback (most recent call last):
  File "hyper_preferential_attachment.py", line 8, in <module>
    import numpy as np
ModuleNotFoundError: No module named 'numpy'
=#  else
        HyperPA.external()
    end
    @test true
end

@test FBDCompare.hypergraphsize((Dict{Int64, Dict}(4 => Dict{Vector{Int64}, Int64}(), 2 => Dict([2, 4] => 1, [2, 6] => 1, [7, 7] => 1), 3 => Dict([5, 6, 6] => 1, [6, 6, 6] => 1), 1 => Dict([3] => 1, [7] => 1)), 1:7)) == 14

FBDCompare.save_figures()


text:
type IIDRandomSequence{S<:ValueSupport,T<:UnivariateDistribution} <: AbstractIDSequence{S}
  d::T
  length::Int64
end
function IIDRandomSequence{T<:UnivariateDistribution}(d::T,length::Int64)
    S = supertype(T).parameters[2]
    return IIDRandomSequence{S,T}(d,length)
end

length(X::IIDRandomSequence) = X.length
distributions(X::IIDRandomSequence) = fill(X.d,length(X))

function insupport(X::IIDRandomSequence,x::AbstractVector)
  length(x)==length(X) && all(insupport(X.d,x))
end

_rand!{T<:Real}(X::IIDRandomSequence, x::AbstractVector{T}) = rand!(X.d,x)

function _logpdf{T<:Real}(X::IIDRandomSequence, x::AbstractVector{T})
  return sum(x_i->logpdf(X.d,x_i),x)
end

mean(X::IIDRandomSequence) = mean(X.d)*ones(length(X))
var(X::IIDRandomSequence) = var(X.d)*ones(length(X))
cov(X::IIDRandomSequence) = Diagonal(var(X))
entropy(X::IIDRandomSequence) = entropy(X.d)*length(X)

function mgf(X::IIDRandomSequence, t::AbstractVector)
  prod([mgf(X.d,t_i) for t_i in t])
end
function cf(X::IIDRandomSequence, t::AbstractVector)
  prod([cf(X.d,t_i) for t_i in t])
end


text:
module Patterns

using Metatheory: binarize, cleanast, alwaystrue
using AutoHashEquals
using TermInterface

"""
Abstract type representing a pattern used in all the various pattern matching backends.
"""
abstract type AbstractPat end

struct UnsupportedPatternException <: Exception
  p::AbstractPat
end

Base.showerror(io::IO, e::UnsupportedPatternException) = print(io, "Pattern ", e.p, " is unsupported in this context")

Base.isequal(a::AbstractPat, b::AbstractPat) = false
TermInterface.arity(p::AbstractPat) = 0
"""
A ground pattern contains no pattern variables and
only literal values to match.
"""
isground(p::AbstractPat) = false
isground(x) = true # literals

# PatVar is equivalent to SymbolicUtils's Slot
"""
    PatVar{P}(name, debrujin_index, predicate::P)

Pattern variables will first match on one subterm
and instantiate the substitution to that subterm.

Matcher pattern may contain pattern variables with attached predicates,
where `predicate` is a function that takes a matched expression and returns a
boolean value. Such a slot will be considered a match only if `f` returns true.

`predicate` can also be a `Type{<:t}`, this predicate is called a
type assertion. Type assertions on a `PatVar`, will match if and only if
the type of the matched term for the pattern variable is a subtype of `T`.
"""
mutable struct PatVar{P} <: AbstractPat
  name::Symbol
  idx::Int
  predicate::P
end
function Base.isequal(a::PatVar, b::PatVar)
  # (a.name == b.name)
  a.idx == b.idx
end
PatVar(var) = PatVar(var, -1, alwaystrue)
PatVar(var, i) = PatVar(var, i, alwaystrue)

"""
If you want to match a variable number of subexpressions at once, you will need
a **segment pattern**.
A segment pattern represents a vector of subexpressions matched.
You can attach a predicate `g` to a segment variable. In the case of segment variables `g` gets a vector of 0 or more
expressions and must return a boolean value.
"""
mutable struct PatSegment{P} <: AbstractPat
  name::Symbol
  idx::Int
  predicate::P
end

PatSegment(v) = PatSegment(v, -1, alwaystrue)
PatSegment(v, i) = PatSegment(v, i, alwaystrues)

"""
Term patterns will match
on terms of the same `arity` and with the same
function symbol `operation` and expression head `exprhead`.
"""
struct PatTerm <: AbstractPat
  exprhead::Any
  operation::Any
  args::Vector
  mod::Module # useful to match against function head symbols and function objs at the same time
  PatTerm(eh, op, args, mod) = new(eh, op, args, mod) #Ref{UInt}(0))
end
TermInterface.istree(::PatTerm) = true
TermInterface.exprhead(e::PatTerm) = e.exprhead
TermInterface.operation(p::PatTerm) = p.operation
TermInterface.arguments(p::PatTerm) = p.args
TermInterface.arity(p::PatTerm) = length(arguments(p))
TermInterface.metadata(p::PatTerm) = p.mod

function TermInterface.similarterm(x::PatTerm, head, args, symtype = nothing; metadata = @__MODULE__, exprhead = :call)
  PatTerm(exprhead, head, args, metadata)
end

# function Base.hash(t::PatTerm, salt::UInt)
#     !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
#     h = t.hash[]
#     !iszero(h) && return h
#     h′ = hash(t.exprhead, hash(t.operation, hash(t.args, salt)))
#     t.hash[] = h′
#     return h′
# end

isground(p::PatTerm) = all(isground, p.args)

# ==============================================
# ================== PATTERN VARIABLES =========
# ==============================================

"""
Collects pattern variables appearing in a pattern into a vector of symbols
"""
patvars(p::PatVar, s) = push!(s, p.name)
patvars(p::PatSegment, s) = push!(s, p.name)
patvars(p::PatTerm, s) = (patvars(operation(p), s); foreach(x -> patvars(x, s), arguments(p)); s)
patvars(x, s) = s

patvars(p) = unique!(patvars(p, Symbol[]))

# ==============================================
# ================== DEBRUJIN INDEXING =========
# ==============================================

function setdebrujin!(p::Union{PatVar,PatSegment}, pvars)
  p.idx = findfirst((==)(p.name), pvars)
end

# literal case
setdebrujin!(p, pvars) = nothing

function setdebrujin!(p::PatTerm, pvars)
  setdebrujin!(operation(p), pvars)
  foreach(x -> setdebrujin!(x, pvars), p.args)
end

#TODO ADD ORIGINAL CODE OF PREDICATE TO PATVAR ?
function to_expr(x::PatVar)
  if x.predicate == alwaystrue
    Expr(:call, :~, x.name)
  else
    Expr(:call, :~, Expr(:(::), x.name, x.predicate))
  end
end

to_expr(x::Any) = x

function to_expr(x::PatSegment)
  Expr(:..., x.predicate == alwaystrue ? Expr(:call, :~, x.name) : Expr(:call, :~, Expr(:(::), x.name, x.predicate)))
end

to_expr(x::PatSegment{typeof(alwaystrue)}) = Expr(:..., Expr(:call, :~, x.name))

to_expr(x::PatSegment{T}) where {T<:Function} = Expr(:..., Expr(:call, :~, Expr(:(::), x.name, nameof(T))))

to_expr(x::PatSegment{<:Type{T}}) where {T} = Expr(:..., Expr(:call, :~, Expr(:(::), x.name, T)))

function to_expr(x::PatTerm)
  pl = operation(x)
  similarterm(Expr(:call, :x), pl, map(to_expr, arguments(x)); exprhead = exprhead(x))
end

# include("rules/patterns.jl")
export AbstractPat
export PatVar
export PatTerm
export PatSegment
export patvars
export setdebrujin!
export isground
export UnsupportedPatternException

end


text:
# TODO I think here is some potential for faster computation
# However, up to 25 states, the algorithm works fine and fast (main knobs are rtol and maxiter)
# This is the specialized version assuming a mass matrix form / linear in dx
# M(x, p)*dx = f(x, p)

# TODO preallocation
function ISInDy(X::AbstractArray, Ẋ::AbstractArray, Ψ::Basis; maxiter::Int64 = 10, rtol::Float64 = 0.99, p::AbstractArray = [], opt::T = ADM()) where T <: DataDrivenDiffEq.Optimise.AbstractSubspaceOptimiser
    @assert size(X)[end] == size(Ẋ)[end]
    nb = length(Ψ.basis)

    # Compute the library and the corresponding nullspace
    θ = Ψ(X, p = p)
    # Init for sweep over the differential variables
    eqs = Operation[]

    @inbounds for i in 1:size(Ẋ, 1)
        dθ = hcat(map((dxi, ti)->dxi.*ti, Ẋ[i, :], eachcol(θ))...)
        Θ = vcat(dθ, θ)
        N = nullspace(Θ', rtol = rtol)
        Q = deepcopy(N) # Deepcopy for inplace

        # Find sparse vectors in nullspace
        # Calls effectively the ADM algorithm with varying initial conditions
        DataDrivenDiffEq.fit!(Q, N', opt, maxiter = maxiter)

        # Compute pareto front
        pareto = map(q->norm([norm(q, 0) ;norm(Θ'*q, 2)], 2), eachcol(Q))
        score, posmin = findmin(pareto)
        # Get the corresponding eqs
        q_best = Q[:, posmin]
        # Remove small entries
        q_best[abs.(q_best) .< opt.λ] .= zero(eltype(q_best))
        rmul!(q_best ,one(eltype(q_best))/maximum(abs.(q_best)))

        # Numerator and Denominator
        # Maybe there is a better way of doing this
        Fn = simplify_constants(simplified_matvec(q_best[nb+1:end], Ψ.basis))
        Fd = simplify_constants(simplified_matvec(q_best[1:nb], Ψ.basis))

        push!(eqs, -Fn/Fd)
    end

    return Basis(eqs, variables(Ψ))
end


text:
using Test

using Tullio: storage_type, promote_storage
using ForwardDiff, FillArrays

@testset "storage_type" begin

    @test storage_type(rand(2), rand(2,3)) == Array{Float64,N} where N
    @test storage_type(rand(2), rand(Float32, 2)) == Vector{Float64}
    @test storage_type(rand(2), rand(Float32, 2,2)) == Array{Float64,N} where N

    Base.promote_type(Matrix{Int}, Vector{Int}) == Array{Int64,N} where N
    Base.promote_type(Matrix{Int}, Matrix{Int32}) == Matrix{Int64}
    Base.promote_type(Matrix{Int}, Vector{Int32}) == Array # != Array{Int64,N} where N
    promote_storage(Matrix{Int}, Vector{Int32}) == Array{Int64,N} where N

    @test storage_type(rand(2), 1:3) == Vector{Float64}
    @test storage_type(rand(Int,2), 1:3) == Vector{Int}
    @test storage_type(1:3.0, 1:3) <: AbstractRange{Float64}

    @test storage_type(rand(2), fill(ForwardDiff.Dual(1,0),2)) == Vector{ForwardDiff.Dual{Nothing,Float64,1}}
    @test storage_type(rand(2), fill(ForwardDiff.Dual(1,0),2,3)) == Array{ForwardDiff.Dual{Nothing,Float64,1}}
    @test storage_type(rand(2), FillArrays.Fill(1.0, 2,2)) == AbstractArray{Float64,N} where N
    @test storage_type(rand(2), FillArrays.Fill(true, 2,2)) == AbstractArray

end

using Tullio: range_expr_walk, divrange, minusrange, subranges, addranges

@testset "range_expr_walk" begin

    for r in [Base.OneTo(10), 0:10, 0:11, 0:12, -1:13]
        for (f, ex) in [
            # +
            (i -> i+1, :(i+1)),
            (i -> i+2, :(i+2)),
            (i -> 3+i, :(3+i)),
            # -
            (i -> -i, :(-i)),
            (i -> i-1, :(i-1)),
            (i -> 1-i, :(1-i)),
            (i -> 2-i, :(2-i)),
            (i -> 1+(-i), :(1+(-i))),
            (i -> -i+1, :(-i+1)),
            (i -> -i-1, :(-i-1)),
            (i -> 1-(2-i), :(1-(2-i))),
            (i -> 1-(-i+2), :(1-(-i+2))),
            # *
            (i -> 2i, :(2i)),
            (i -> 2i+1, :(2i+1)),
            (i -> -1+2i, :(-1+2i)),
            (i -> 1-3i, :(1-3i)),
            (i -> 1-3(i+4), :(1-3(i+4))),
            # ÷
            (i -> i÷2, :(i÷2)),
            (i -> 1+i÷3, :(1+i÷3)),
            # triple...
            (i -> i+1+2, :(i+1+2)),
            (i -> 1+2+i, :(1+2+i)),
            (i -> 2i+3+4, :(2i+3+4)),
            (i -> 1+2+3i+4, :(1+2+3i+4)),
            (i -> 1+2+3+4(-i), :(1+2+3+4(-i))),
            # evil
            (i -> (2i+1)*3+4, :((2i+1)*3+4)),
            (i -> 3-(-i)÷2, :(3-(-i)÷2)), # needs divrange_minus
            ]
            rex, i = range_expr_walk(:($r .+ 0), ex)
            @test issubset(sort(f.(eval(rex))), r)
        end
        @test minusrange(r) == divrange(r, -1)

        @test issubset(subranges(r, 1:3) .+ 1, r)
        @test issubset(subranges(r, 1:3) .+ 3, r)
        @test union(subranges(r, 1:3) .+ 1, subranges(r, 1:3) .+ 3) == r

        @test issubset(addranges(r, 1:3) .- 1, r)
        @test issubset(addranges(r, 1:3) .- 3, r)
        @test sort(union(addranges(r, 1:3) .- 1, addranges(r, 1:3) .- 3)) == r
    end
end

using Tullio: cleave, trisect, productlength

@testset "threading" begin
    @test cleave((1:10, 1:4, 7:8)) == ((1:4, 1:4, 7:8), (5:10, 1:4, 7:8))
    @test cleave((7:8, 9:9)) == ((7:7, 9:9), (8:8, 9:9))
    @test cleave((1:4,)) == ((1:2,), (3:4,))
    @test cleave(()) == ((), ())

    @test trisect((1:9, 11:12)) == ((1:3, 11:12), (4:6, 11:12), (7:9, 11:12))
    @test trisect((1:9,)) == ((1:3,), (4:6,), (7:9,))
    @test trisect(()) == ((), (), ())

    @test sum(productlength, trisect((1:10, 11:20))) == 100

    for r1 in [1:10, 3:4, 5:5], r2 in [11:21, -3:-2, 0:0], r3 in [2:7, 8:9, 0:0]
        tup = (r1,r2,r3)
        len = productlength(tup)
        @test len == sum(productlength, cleave(tup))
        @test len == sum(productlength, trisect(tup))
    end
end

using Tullio: @capture_

@testset "capture_ macro" begin
    EXS  = [:(A[i,j,k]),  :(B{i,2,:}),  :(C.dee), :(fun(5)),   :(g := h+i),        :(k[3] += l[4]), :([m,n,0]) ]
    PATS = [:(A_[ijk__]), :(B_{ind__}), :(C_.d_), :(f_(arg_)), :(left_ := right_), :(a_ += b_),     :([emm__]) ]
    # @test length(EXS) == length(PATS)
    @testset "ex = $(EXS[i])" for i in eachindex(EXS)
        for j in eachindex(PATS)
        @eval res = @capture_($EXS[$i], $(PATS[j]))
        if i != j
            @test res == false
        else
            @test res == true
            if i==1
                @test A == :A
                @test ijk == [:i, :j, :k]
            elseif i==3
                @test C == :C
                @test d == :dee
            elseif i==5
                @test left == :g
                @test right == :(h+i)
            elseif i==7
                @test emm == [:m, :n, 0]
            end
        end
        end
    end

end

using Tullio: leibnitz
using ForwardDiff

@testset "symbolic gradients" begin

    @testset "ex = $ex" for ex in [
        :(x*y + 1/z),
        :(x*y + z*x*z^2),
        :((x/y)^2 + inv(z)),
        :((x+2y)^z),
        :(1/(x+1) + 33/(22y) -4/(z/4)),
        :(inv(x+y) + z^(-2)),

        :(sqrt(x) + 1/sqrt(y+2z)),
        :(inv(sqrt(x)*sqrt(y)) + sqrt(2*inv(z))),
        :(x * z / sqrt(y * z)),

        :(log(x/y) - log(z+2)),
        :(log(x*y*z) - 33y),

        :(2exp(x*y*z)),
        :(exp((x-y)^2/2)/z),
    ]

        dfdx = leibnitz(ex, :x)
        dfdy = leibnitz(ex, :y)
        dfdz = leibnitz(ex, :z)

        @eval f_(x,y,z) = $ex
        @eval f_x(x,y,z) = $dfdx
        @eval f_y(x,y,z) = $dfdy
        @eval f_z(x,y,z) = $dfdz

        xyz = rand(Float32, 3)

        # check correctness
        gx, gy, gz = ForwardDiff.gradient(xyz -> f_(xyz...), xyz)
        @test f_x(xyz...) ≈ gx
        @test f_y(xyz...) ≈ gy
        @test f_z(xyz...) ≈ gz

        # don't accidentally make Float64
        @test 0f0 + f_x(xyz...) isa Float32
        @test 0f0 + f_y(xyz...) isa Float32
        @test 0f0 + f_z(xyz...) isa Float32

    end

end

macro cse(ex)
    esc(Tullio.commonsubex(ex))
end

@testset "common subexpressionism" begin

    x,y,z = 1.2, 3.4, 5.6

    @test (x+y)*z/(x+y) ≈
        @cse (x+y)*z/(x+y)

    @test x*y*z + 2*x*y*z ≈
        @cse x*y*z + 2*x*y*z

    @test (sqrt(inv(x)) * inv(sqrt(y)) + inv(x)/inv(z)) ≈
        @cse (sqrt(inv(x)) * inv(sqrt(y)) + inv(x)/inv(z))

    @test (a1 = inv(x); b1 = inv(x); c1 = inv(x)*inv(y)) ≈
        @cse (a = inv(x); b = inv(x); c = inv(x)*inv(y))

    # setting a, b (outside @test)
    (a1 = inv(x); b1 = inv(x); c1 = inv(x)*inv(y))
    @cse (a = inv(x); b = inv(x); c = inv(x)*inv(y))

    @test a1 ≈ a
    @test b1 ≈ b
    @test c1 ≈ c

    # updating a, b, etc
    a = 1
    @cse a = a + x*y/z
    @test a ≈ 1 + x*y/z

    a = 1
    b = 1
    @cse begin
        a = a + x*y/z
        b += x*y + z
    end
    @test a ≈ 1 + x*y/z
    @test b ≈ 1 + x*y + z

    a = a1 = 1
    b = b1 = 1
    begin
        θ = inv(x)*inv(y+a)
        a1 += θ^2
        b1 = θ^3 + inv(x)
    end
    @cse begin
        θ = inv(x)*inv(y+a)
        a += θ^2
        b = θ^3 + inv(x)
    end
    @test a ≈ a1
    @test b ≈ b1

end


text:
# Julia wrapper for header: /usr/include/libavutil/avutil.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    #avutil_version,
    avutil_configuration,
    avutil_license,
    av_get_picture_type_char

function avutil_version()
    ccall((:avutil_version,libavutil),Uint32,())
end

function avutil_configuration()
    ccall((:avutil_configuration,libavutil),Ptr{Uint8},())
end

function avutil_license()
    ccall((:avutil_license,libavutil),Ptr{Uint8},())
end

function av_get_picture_type_char(pict_type::AVPictureType)
    ccall((:av_get_picture_type_char,libavutil),Uint8,(AVPictureType,),pict_type)
end


text:
const tolerance = 10.0^-10

struct Line
    x1::Number
    y1::Number

    x2::Number
    y2::Number
end

Base.:(==)(lhs::Line, rhs::Line) = lhs.x1 == rhs.x1 && lhs.y1 == rhs.y1 && lhs.x2 == rhs.x2 && lhs.y2 == rhs.y2

# How much is missing to get to the closest new "grid square"
function rayDelta(n, a)
    s = sign(a)

    if s > 0
        return floor(n + 1) - n

    elseif s < 0
        return ceil(n - 1) - n

    else
        return 0
    end
end

function pointsOnLine(line::Line)
    x1, y1, x2, y2 = line.x1, line.y1, line.x2, line.y2
    res = Set{Tuple{Number, Number}}()
    theta = atan(y2 - y1, x2 - x1)

    x, y = x1, y1

    velX = cos(theta)
    velY = sin(theta)

    velX = isapprox(velX, 0, atol=tolerance) ? 0 : velX
    velY = isapprox(velY, 0, atol=tolerance) ? 0 : velY

    while velX != 0 && !isapprox(x, x2, atol=tolerance) || velY != 0 && !isapprox(y, y2, atol=tolerance)
        dx = rayDelta(x, velX)
        dy = rayDelta(y, velY)

        a = isnan(dx / velX) ? Inf : dx / velX
        b = isnan(dy / velY) ? Inf : dy / velY

        if a < b
            x += dx
            y += a * velY

        else
            x += b * velX
            y += dy
        end

        fx = floor(Int, x)
        fy = floor(Int, y)

        push!(res, (fx, fy))
    end

    push!(res, (x1, y1))
    push!(res, (x2, y2))

    return collect(res)
end