text:
module PNG16x16

using Colors
using Base64

export write_png_as_data, write_png

function write_png_as_data(io::IO, cs::AbstractMatrix{<:Color})
    write(io, "data:image/png;base64,")
    b64enc = Base64EncodePipe(io)
    write_png(b64enc, cs)
    close(b64enc)
end

function write_png(io::IO, cs::AbstractMatrix{<:Color})
    buf = IOBuffer() # to calculate chunk CRCs
    n = 16 # 16 x 16
    u8(x) = write(buf, UInt8(x & 0xFF))
    u16(x) = (u8((x & 0xFFFF)>>8); u8(x))
    u32(x) = (u16((x & 0xFFFFFFFF)>>16); u16(x))
    b(bstr) = write(buf, bstr)
    function palette(c::Color)
        rgb24 = convert(RGB24,c)
        u8(rgb24.color>>16); u8(rgb24.color>>8); u8(rgb24.color)
    end
    crct(x) = (for i = 1:8; x = x & 1==1 ? 0xEDB88320 ⊻ (x>>1) : x>>1 end; x)
    table = UInt32[crct(i) for i = 0x00:0xFF]
    function crc32()
        seekstart(buf)
        crc = 0xFFFFFFFF
        while !eof(buf)
            crc = (crc>>8) ⊻ table[(crc&0xFF) ⊻ read(buf, UInt8) + 1]
        end
        u32(crc ⊻ 0xFFFFFFFF)
    end
    flush() = write(io, take!(seekstart(buf)))

    # The following is a pre-encoded 256-indexed-color PNG with size of 16x16.
    # We only rewrite "pallets".
    b(b"\x89PNG\x0D\x0A\x1A\x0A")
    # Image header
    u32(13); flush(); b(b"IHDR"); u32(n); u32(n); u8(8); u8(3); u8(0); u8(0); u8(0); crc32()
    # Palette
    u32(n * n * 3); flush();
    b(b"PLTE")
    for y = 1:n, x = 1:n
        palette(cs[y,x])
    end
    crc32()
    # Image data
    u32(58); flush(); b(b"IDAT")
    b(b"\x78\xDA\x63\x64\x60\x44\x03\x02\xE8\x02\x0A\xE8\x02\x06\xE8\x02")
    b(b"\x0E\xE8\x02\x01\xE8\x02\x09\xE8\x02\x05\xE8\x02\x0D\xE8\x02\x13")
    b(b"\xD0\x05\x16\xA0\x0B\x6C\x40\x17\x38\x80\x2E\x70\x01\x5D\xE0\x01")
    b(b"\xBA\xC0\x07\x34\x3E\x00\x54\x4D\x08\x81"); crc32()
    # Image trailer
    u32(0); flush(); b(b"IEND"); crc32()
    flush()
end

"""
# Image data
using CodecZlib
raw = IOBuffer()
for y = 0:15
    write(raw, UInt8(1)) # filter: SUB
    write(raw, UInt8(y*16)) # line head
    write(raw, UInt8[1 for i=1:15]) # left + 1
end
flush(raw)
cd = ZlibCompressorStream(raw,level=9)
flush(cd)
seekstart(cd)
@show read(cd) # UInt8[0x78, 0xda, 0x63, 0x64, ...
"""

end # module


text:
include("header.jl")

#=
# Velocity profiles
In this notebook, we will discuss a useful way to visualize a flow in fluid
mechanics: the **velocity profile**.
=#
#-

# ### Set up the module
using MAE103
#-
using Plots

#=
To start, let's consider one of the most basic types of flow: the linear shear flow,
often called **Couette flow**. This type of flow is generated, for example, between
two parallel walls when one wall is moving and the other is stationary. In the figure
below, the top wall is moving with velocity $U$, the lower wall is stationary:

Because of the **no-slip condition**, the fluid next to each wall moves with it. The
fluid next to the upper wall moves at velocity $U$, the fluid next to the lower wall
is at rest.

The velocity profile depicts $u(y)$, the horizontal component of velocity as a
function of vertical position, $y$. This increases linearly from the lower to the
upper wall. In fact, the function is just

$$u(y) = Uy/H$$

The arrows indicate the direction that the fluid is moving, and the lengths of
the arrows indicate the relative speed at that $y$ position.
=#
#-
#=
### Shear stress
We can calculate the viscous shear stress at any location with the equation

$$\tau = \mu \frac{\mathrm{d}u}{\mathrm{d}y}$$

In fact, in Couette flow, shear stress is the same everywhere:

$$\tau = \mu U/H$$

This represents, e.g., the resistive force (per unit area) that the fluid applies
(to the left) on the upper wall while the wall moves to the right.

It also represents the force per unit area that the lower wall must withstand in
order to remain stationary. The fluid is trying to drag the lower wall to the right,
so this must be balanced by some reaction force in the wall to the left.

The **rate of work** (per unit area) required to move the upper wall is $\dot{W} = \tau U$.
This would be measured in units of, e.g., $W/m^2$.
=#
#-
#=
### Plotting velocity profiles
The arrows are helpful, but you can also plot a velocity profile without them.
For example, consider the following velocity:

$$u(y) = \frac{4U_c}{H^2} y (H - y)$$

The coefficient $U_c$ is a speed, and $H$ is the gap height.
=#
#=
Let's define a function that evaluates this velocity. Here, `y`, `Uc`, and `H`
are to be given as arguments to the function.
=#
u(y,Uc,H) = 4*Uc/H^2*y*(H-y)

#=
Suppose the gap height $H$ is 1 cm and the speed $U_c$ is 1 m/s. We will evaluate
this velocity at a range of locations between 0 and $H$. Here, we evaluate
the velocity at 101 points. There's nothing magical about this number; it just
is enough to make the curve look smooth. Also, note that we have to give
units for the `y` values, as well:
=#
H = 1u"cm"
Uc = 1u"m/s"
y = range(0u"cm",H,length=101)

#=
Now we evaluate the velocity function at the range of $y$ locations. (Remember that
the `.` vectorizes the evaluation of a function. Here, we are telling it to
evaluate the function `u` at each element of the `y` array.) The resulting
velocities are assigned to `v`.
=#
v = u.(y,Uc,H)

#=
Notice that $u$ is 0 at the beginning and end of the range. Let's plot it. But
let's plot it as a velocity profile, which means we make velocity the 'x' axis and $y$
the 'y' axis:
=#
plot(v,y,xlim=(0,Inf),ylim=(0,Inf),xlabel=L"u(y)",ylabel=L"y")

#=
NOTE: The use of `xlim=(0,Inf),ylim=(0,Inf)` simply ensures that the axes start at 0
and end at the upper range of the data in `v` and `y`, respectively.

The top and bottom of this plot suggest that these are stationary walls where the
flow is at rest. In fact, this is the velocity profile associated with pressure-driven
flow through the gap.
=#
#-
#=
Suppose we wish to plot the velocity in different units. We can easily do this with
the `xunit` flag. For example, in `ft/s`:
=#
plot(v,y,xunit=u"ft/s",xlim=(0,Inf),ylim=(0,Inf),xlabel=L"u(y)",ylabel=L"y")
#-
#=
In the previous notebook, we showed that we can set velocity using
the `Velocity` function. We can do that here if we like. We'll call it `vv`.
We have to use the `.` notation to make sure each element in `v` gets the
same treatment:
=#
vv = Velocity.(v)
plot(vv,y,xlim=(0,Inf),ylim=(0,Inf),xlabel="u(y)",ylabel="y")

#=
The `Velocity` function simply *wraps* the velocity data in `v`. It is
probably not immediately clear why we would need to do this, since the
data are already in the correct units (m/s). However, it
will be useful later when we use functions that perform calculations,
and we wish to make sure we don't try to provide the wrong type of data
to those functions.
=#
#-

#=
#### Viscometry
We can use the relationship between the velocity and shear stress for many practical
tasks. In particular, it can be used to measure the viscosity of a given fluid,
in a device called a **viscometer**. This device consists of a rotating inner
cylinder of radius $R$ with prescribed angular velocity $\Omega$ and a stationary
outer wall, and a fluid with unknown viscosity, $\mu$, in the gap of size $H$ between
them.

The basic idea of a viscometer is that one knows the angular velocity $\Omega$
and measures the *torque*, $T$, required to rotate it. This torque is directly
attributable to the viscous shear stress acting on the inner cylinder.

Even though this is an annular region and not a straight gap like in the basic
Couette flow, the curvature does not matter here because the gap $H$ is very small
compared to the radius $R$.

Comparing this gap flow with the Couette flow above, the wall velocity here is
$U = \Omega R$, so the wall shear stress applied on the inner cylinder is

$$\tau = \mu \Omega R/H$$

How do we relate this to the measured torque? Torque is the force applied on the
wall times the distance to the center, $R$. The shear stress is uniformly applied
on the inner cylinder, so the force applied by the fluid is the shear stress times
the surface area,

$$T = \tau (2\pi R L) R$$

The term in parentheses is the surface area of the inner cylinder, where $L$ is its
length (into the page).

Thus, putting all of this together:

$$T = 2\pi L\mu \Omega R^3/H$$
=#
#-
#=
#### A viscometer example
Let's suppose that we take measurements of a viscometer with the following
characteristics: inner radius $R = 6$ cm, gap height $H = 0.125$ cm, and length
$L = 13$ cm. We will convert everything to SI (i.e., meters):
=#
R = 6u"cm"
H = 0.125u"cm"
L = 13u"cm";

#=
We run the viscometer at 6 different speeds and measure the required torque in
each case. Here is the data we measured. Angular velocity is given in rad/s,
and torque in N$\cdot$m
=#
Ω = [1.0,2.0,3.0,4.0,5.0,6.0]u"rad/s"

T = [0.181,0.36,0.54,0.72,0.899,1.086]u"N*m";

#-
#=
Let's plot the data to see the trend qualitatively:
=#
scatter(Ω,T,xlim=(0,Inf),ylim=(0,Inf),xlabel="Angular velocity",ylabel="Torque")

#=
It looks quite linear! This is expected, since torque is linearly dependent on
the angular velocity in the relationship above.
Based on the equation above, we expect that the ratio of these two quantities
will be equal to $2\pi L\mu R^3/H$. So if we multiply the ratio by $H/(2\pi L R^3)$,
we expect

$$\dfrac{T H}{2\pi L R^3 \Omega} = \mu$$

Let's do this for our data. Since we are doing this for each data point,
we will get an array of values. (Notice the `.` in front of `/` when we divide by `Ω`.
Dividing by an array needs to be done in a vectorized manner.) We assign this
array the name `a`:
=#
a = T./Ω*H/(2π*L*R^3)

#=
We expect `a` to give us viscosity, but there's always the chance for
experimental error. So we will take the mean of `a`. And we will use the
`Viscosity` function to put this in the default units of viscosity:
=#
Viscosity(mean(a))

#=
It's quite viscous! Only slightly smaller viscosity than glycerin.
=#
Viscosity(Glycerin)

ext: jl


text:
#                     5
#                     |
# 6 -> 7 -> 4 -> 8 -> 9
#                     |
#           1 -> 2 -> 3

s1 = EuclideanStation(1, EuclideanCoordinate(3, 1), "test", false)
s2 = EuclideanStation(2, EuclideanCoordinate(4, 1), "test", false)
s3 = EuclideanStation(3, EuclideanCoordinate(5, 1), "test", false)
s5 = EuclideanStation(5, EuclideanCoordinate(5, 4), "test", false)

s9 = EuclideanStation(9, EuclideanCoordinate(5, 3), "test", false)

s6 = EuclideanStation(6, EuclideanCoordinate(1, 3), "test", false)
s7 = EuclideanStation(7, EuclideanCoordinate(2, 3), "test", false)
s4 = EuclideanStation(4, EuclideanCoordinate(3, 3), "test", false)
s8 = EuclideanStation(8, EuclideanCoordinate(4, 3), "test", false)

line1 = Line("S1")
line2 = Line("S2")
e1 = ProcessedEdge(s1, s2, line1, 0, 1, false)
e2 = ProcessedEdge(s2, s3, line1, 0, 1, false)
e3 = ProcessedEdge(s3, s9, line1, 2, 1, false)
e4 = ProcessedEdge(s9, s5, line1, 2, 1, false)

e5 = ProcessedEdge(s6, s7, line2, 0, 1, false)
e6 = ProcessedEdge(s7, s4, line2, 0, 1, false)
e7 = ProcessedEdge(s4, s8, line2, 0, 1, false)
e8 = ProcessedEdge(s8, s9, line2, 0, 1, false)
e9 = ProcessedEdge(s9, s5, line2, 0, 1, false)
transit_map = InputGraph([s1, s2, s3, s4, s5,
                          s6, s7, s8, s9],
    [e1, e2, e3, e4,
     e5, e6, e7, e8, e9],
    [line1, line2])

reduced_map = reduce_transitmap(transit_map)

@test length(reduced_map.stations) == 4
@test length(reduced_map.edges) == 4

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
# ------------ #
# Main routine #
# ------------ #

const DEFAULT_SOLVER = IpoptSolver(hessian_approximation="limited-memory",
                                   print_level=2)

"""
TODO: Document the rest of the arguments

`mf` should be a function that computes the empirical moments of the
model. It can have one of two call signatures:

1. `mf(θ)`: computes moments, given only a parameter vector
2. `mf(θ, data)`: computes moments, given a parameter vector and an
   arbitrary object that contains the data necessary to compute the
   moments. Examples of `data` is a matrix, a Dict, or a DataFrame.
   The data argument is not used internally by these routines, but is
   simply here for user convenience.

The `mf` function should return an object of type Array{Float64, 2}
"""
function gmm(mf::Function, theta::Vector, W::Array{Float64, 2},
             instruments::Union{Void, Matrix{Float64}}=nothing;
             solver=DEFAULT_SOLVER, data=nothing,
             mgr::IterationManager=OneStepGMM())
    npar = length(theta)
    theta_l = fill(-Inf, npar)
    theta_u = fill(+Inf, npar)
    gmm(mf, theta, theta_l, theta_u, W,  instruments, solver=solver, data=data,
        mgr=mgr)
end

function gmm(mf::Function, theta::Vector, theta_l::Vector, theta_u::Vector,
             W::Array{Float64, 2},
             instruments::Union{Void, Matrix{Float64}}=nothing;
             solver=DEFAULT_SOLVER,
             data=nothing,
             mgr::IterationManager=OneStepGMM())

    # NOTE: all handling of data happens right here, because we will use _mf
    #       internally from now on.
    if max_args(mf) == 1
        _mf(theta) = mf(theta)
    else
        _mf(theta) = mf(theta, data)
    end

    # Instrument handling here... use _mfi internally afterwards
    if instruments === nothing
        _mfi(theta) = _mf(theta)
    else
        _mfi(theta) = row_kron(_mf(theta), instruments)
    end

    mf0        = _mfi(theta)
    nobs, nmom = size(mf0)
    npar       = length(theta)

    nl         = length(theta_l)
    nu         = length(theta_u)

    @assert nl == nu
    @assert npar == nl
    @assert nobs > nmom
    @assert nmom >= npar

    ## mf is n x m
    smf(theta) = reshape(sum(_mfi(theta),1), nmom, 1);
    _smf(theta) = vec(smf(theta))
    function smf!(gg, θ::AbstractVector)
        gg[:] = smf(θ)
        gg
    end

    Dsmf = ForwardDiff.jacobian(_smf, mutates=false)

    l = theta_l
    u = theta_u
    lb = Float64[]
    ub = Float64[]

    # begin iterations
    ist = IterationState(0, 10.0, theta)

    # Define these outside while loop so they are available after it
    NLPE = GMMNLPE(_mfi, smf, Dsmf, mgr, W)
    m = SolverInterface.NonlinearModel(solver)

    while !(finished(mgr, ist))
        NLPE = GMMNLPE(_mfi, smf, Dsmf, mgr, W)
        m = SolverInterface.NonlinearModel(solver)
        SolverInterface.loadproblem!(m, npar, 0, l, u, lb, ub, :Min, NLPE)
        SolverInterface.setwarmstart!(m, theta)
        SolverInterface.optimize!(m)

        # update theta and W
        theta = SolverInterface.getsolution(m)
        W = optimal_W(_mfi, theta, mgr.k)

        # update iteration state
        ist.n += 1
        ist.change = maxabs(ist.prev - theta)
        ist.prev = theta
    end

    r = GMMResult(SolverInterface.status(m),
                  SolverInterface.getobjval(m),
                  SolverInterface.getsolution(m),
                  nmom, npar, nobs)
    GMMEstimator(NLPE, r)
end


text:
using Pkg
using GLMakie, Test
path = normpath(joinpath(dirname(pathof(AbstractPlotting)), "..", "test", "ReferenceTests"))
Pkg.develop(PackageSpec(path = path))
using ReferenceTests
files, recorded = ReferenceTests.record_tests()
recorded = ReferenceTests.basedir("recorded")
ReferenceTests.reference_tests(recorded)
# needs GITHUB_TOKEN to be defined
# ReferenceTests.upload_reference_images()
# Run the below, to generate a html to view all differences:
# recorded, ref_images, scores = ReferenceTests.reference_tests(recorded)
# ReferenceTests.generate_test_summary("preview.html", recorded, ref_images, scores)


text:
module Stars
# using Stars

using GDAL
using ArchGDAL; const AG = ArchGDAL

using CoordinateTransformations
export AffineMap;

using StaticArrays
using GeoFormatTypes
const GFT = GeoFormatTypes
# include("shp/GeoDataFrames.jl")

import Statistics
import Statistics: mean

include("tools_Ipaper.jl")
include("geo_basic.jl")
include("geoarray.jl")

include("st_bbox.jl")
include("st_crs.jl")

include("st_dim.jl")
include("st_coordinates.jl")

include("tools_ratser.jl")

include("utils/utils.jl")
# include("raster/Raster.jl")
include("raster/ncread.jl")
include("raster/clip.jl")

include("gdal_basic.jl")
# include("gdal_info.jl")
include("gdal_polygonize.jl")

include("st_plot.jl")
include("st_subset.jl")

include("st_read.jl")

include("st_write.jl")
include("st_as_sf.jl")
include("st_as_stars.jl")
include("sp/sp.jl")

include("st_mask.jl")
# include("hydro/hydro.jl")
# include("hydro/snap_pour_points.jl")

end


text:
using GPUifyLoops
using Test

kernel(A::Array) = kernel(Val(:CPU), A)
function kernel(::Val{Dev}, A) where Dev
    @setup Dev

    @loop for i in (1:size(A,1);
                    threadIdx().x)
        A[i] = 2*A[i]
    end
    @synchronize
end

@testset "Array" begin
    data = Array{Float32}(undef, 1024)
    kernel(data)
end

@static if Base.find_package("CuArrays") !== nothing
    using CuArrays
    using CUDAnative

    @eval function kernel(A::CuArray)
        @cuda threads=length(A) kernel(Val(:GPU), A)
    end

    @testset "CuArray" begin
        data = CuArray{Float32}(undef, 1024)
        kernel(data)
    end
end

# Scratch arrays

function f1()
    @setup :GPU
    A = @scratch Int64 (12, 3) 2
    @test A.data isa GPUifyLoops.MArray
    @test size(A.data) == (1,)
end

function f2()
    @setup :GPU
    A = @scratch Int64 (12, 3) 1
    @test A.data isa GPUifyLoops.MArray
end

function f3()
    @setup :CPU
    A = @scratch Int64 (12, 3) 1
    @test A isa GPUifyLoops.MArray
end

@testset "Scratch Arrays" begin
    f1()
    f2()
    f3()
end


text:
using TermWin

s = open( readall, joinpath( Pkg.dir( "TermWin" ), "src", "TermWin.jl" ) )
tshow( s )


text:
import CSV

function table(path::String; args...) :: String
  table(CSV.file(path))
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

num_of_rays = 0
small_image = HDRimage(1, 1)
camera = OrthogonalCamera(1.0)
tracer = ImageTracer(small_image, camera, 10, PCG())

function trace_ray(ray::Ray)
    # num_of_rays
    point = at(ray, 1.)

    # Check that all the rays intersect the screen within the region [−1, 1] × [−1, 1]
    @test abs(point.x) < 1e-5
    @test -1.0 <= point.y <= 1.0
    @test -1.0 <= point.z <= 1.0

    global num_of_rays += 1

    return RGB{Float64}(0.0, 0.0, 0.0)
end

fire_all_rays!(tracer, trace_ray)

# Check that the number of rays that were fired is what we expect (10²)
@test num_of_rays == 100


text:
GP_ENV["ALLFIGS"] = Dict{String, Figure}()
GP_ENV["CURFIG"]  = nothing
GP_ENV["CURAXES"] = nothing

"""
    gcf()

Return the current active Figure or a new figure if there isn't one.
See also [`gca`](@ref).
"""
gcf() = isdef(GP_ENV["CURFIG"]) ? GP_ENV["CURFIG"] : Figure() # no ifelse here

"""
    gca()

Return the current active Axes and `nothing` if there isn't one.
See also [`gcf`](@ref).
"""
gca() = GP_ENV["CURAXES"] # if nothing, whatever called it will create

"""
    check_axes(a)

Internal function to check if `a` is defined, if not it calls `gca` if `gca` also
returns `nothing`, it adds axes. Used in `plot!` etc.
"""
function check_axes(a::Option{Axes}; dims=2)
    isdef(a) || (a = gca())
    if isdef(a)
        # check if has the right dims otherwise overwrite
        if dims == 3
            if isa(a, Axes3D)
                return a
            else
                f = parent(a)
                ha = hash(a)
                i = findlast(e->hash(e) === ha, f.axes)
                f.axes[i] = Axes3D{get_backend(f)}(parent=f.id)
            end
        else
            if isa(a, Axes2D)
                return a
            else
                f = parent(a)
                ha = hash(a)
                i = findlast(e->hash(e) === ha, f.axes)
                f.axes[i] = Axes2D{get_backend(f)}(parent=f.id)
            end
        end
    else
        if dims == 3
            return add_axes3d!()
        else
            return add_axes2d!()
        end
    end
end

"""
    get_backend(f)

Return the backend type associated with figure `f`.
"""
get_backend(f::Figure{B}=gcf()) where {B} = B

"""
    palette(v::Vector{Color})

Set the default palette. To see the current one, write `GPlot.GP_ENV["PALETTE"]` in the REPL.
"""
set_palette(v::Vector{<:Color}) = (GP_ENV["PALETTE"] = v; GP_ENV["SIZE_PALETTE"] = length(v); ∅)

"""
    continuous_preview(b)

Set the continuous preview on (`b=true`, default value) or off.
"""
continuous_preview(b::Bool) = (GP_ENV["CONT_PREVIEW"] = b)


text:
using Pkg
Pkg.activate((@__DIR__)*"/..")

using Plasmo
using MPI
using PipsSolver

MPI.Init()
comm = MPI.COMM_WORLD
ncores = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

Ns = 8
SPP = round(Int, floor(Ns/ncores))

function create_simple_model(d)
    node = OptiNode()
    @variable(node, 0<=p<=10, start=5)
    @variable(node, u)
    @variable(node, x)
    @constraint(node, x >= p)
    @constraint(node, p + u == d)
    return node
end

#Setup processor information
demand = rand(Ns)*10
graph = OptiGraph()

first_stage = @optinode(graph)
@variable(first_stage,0 <= x <= 8)
@objective(first_stage,Min,x)

#Add the master model to the graph
# master_node = add_node!(graph,master)
# scenm=Array{JuMP.Model}(undef,Ns)
#scen_nodes = Array{ModelNode}(undef,Ns)
owned = []
s = 1
subgraph = OptiGraph()
add_subgraph!(graph,subgraph)
#split scenarios between processors
for j in 1:Ns
    global s
    if round(Int, floor((s-1)/SPP)) == rank
        push!(owned, s)
        node = create_simple_model(demand[j])
        add_node!(subgraph,node)

        #connect children and parent variables
        @linkconstraint(graph, first_stage[:x] == node[:x])

        #reconstruct second stage objective
        @objective(node,Min,1/Ns*(node[:p] + 3*node[:u]))
    else #Ghost nodes
        node = add_node!(subgraph)
    end
    s = s + 1
end

#create a link constraint between the subproblems (PIPS-NLP supports this kind of constraint)
@linkconstraint(graph, (1/Ns)*sum(getnode(subgraph,s)[:p] for s in owned) == 8)

if rank == 0
    println("Solving with PIPS-NLP")
end

pipsnlp_solve(graph)

MPI.Finalize()


text:
using Documenter, MemorableUniqueIdentifier

makedocs(
    modules=[MemorableUniqueIdentifier],
    sitename="Memorable Unique Identifier"
)

max_stars_repo_name: rustyconover/muid-julia

max_issues_repo_name: rustyconover/muid-julia

max_forks_repo_name: rustyconover/MemorableUniqueIdentifier


text:
using RandomLineChart
using Documenter

DocMeta.setdocmeta!(RandomLineChart, :DocTestSetup, :(using RandomLineChart); recursive=true)

makedocs(;
    modules=[RandomLineChart],
    authors="arjun1638 <arjun1638@gmail.com> and contributors",
    repo="https://github.com/arjun1638/RandomLineChart.jl/blob/{commit}{path}#{line}",
    sitename="RandomLineChart.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://arjun1638.github.io/RandomLineChart.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/arjun1638/RandomLineChart.jl",
    devbranch="main",
)


text:
using Test

using AnalyzeRegistry

@testset "AnalyzeRegistry" begin
    general = general_registry()
    @test isdir(general)
    @test all(isdir, find_packages())
    # Test some properties of the `Measurements` package.  NOTE: they may change
    # in the future!
    measurements = analyze(joinpath(general, "M", "Measurements"))
    @test measurements.reachable
    @test measurements.docs
    @test measurements.runtests
    @test !measurements.buildkite
    # Test results of a couple of packages.  Same caveat as above
    packages = [joinpath(general, p...) for p in (("C", "Cuba"), ("P", "PolynomialRoots"))]
    results = analyze(packages)
    cuba, polyroots = results
    @test length(filter(p -> p.reachable, results)) == 2
    @test length(filter(p -> p.runtests, results)) == 2
    @test cuba.drone
    @test !polyroots.docs # Documentation is in the README!
    # We can also use broadcasting!
    @test Set(results) == Set(analyze.(packages))
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using ForwardDiff

""" Find segment from slave element corresponding to master element nodes.

Parameters
----------
x1_, n1_
    slave element geometry and normal direction
x2
    master element node to project onto slave

Returns
-------
xi
    dimensionless coordinate on slave corresponding to
    projected master

"""
function project_from_master_to_slave{E<:MortarElements2D}(
    slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector;
    tol=1.0e-10, max_iterations=20, debug=false)

    x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
    dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
    n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
    dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
    cross2(a, b) = cross([a; 0], [b; 0])[3]
    R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
    dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))

    xi1 = 0.0
    xi1_next = 0.0
    dxi1 = 0.0
    for i=1:max_iterations
        dxi1 = -R(xi1)/dR(xi1)
        dxi1 = clamp(dxi1, -0.3, 0.3)
        xi1_next = clamp(xi1 + dxi1, -1.0, 1.0)
        if norm(xi1_next - xi1) < tol
            return xi1_next
        end
        if debug
            info("xi1 = $xi1")
            info("R(xi1) = $(R(xi1))")
            info("dR(xi1) = $(dR(xi1))")
            info("dxi1 = $dxi1")
            info("norm = $(norm(xi1_next - xi1))")
            info("xi1_next = $xi1_next")
        end
        xi1 = xi1_next
    end

    info("x1 = $x1_")
    info("n1 = $n1_")
    info("x2 = $x2")
    info("xi1 = $xi1, dxi1 = $dxi1")
    info("-R(xi1) = $(-R(xi1))")
    info("dR(xi1) = $(dR(xi1))")
    error("find projection from master to slave: did not converge")

end

function project_from_slave_to_master{E<:MortarElements2D}(
    master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI;
    tol=1.0e-10, max_iterations=20)

    x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
    dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
    cross2(a, b) = cross([a; 0], [b; 0])[3]
    R(xi2) = cross2(x2(xi2)-x1, n1)
    dR(xi2) = cross2(dx2(xi2), n1)

    xi2 = 0.0
    dxi2 = 0.0
    for i=1:max_iterations
        dxi2 = -R(xi2) / dR(xi2)
        xi2 += dxi2
        if norm(dxi2) < tol
            return xi2
        end
    end

    error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")

end

"""
Frictionless 2d finite sliding contact with forwarddiff.

true/false flags: finite_sliding, friction, use_forwarddiff
"""
function assemble!(problem::Problem{Contact}, time::Float64,
                   ::Type{Val{1}},     ::Type{Val{true}},
                   ::Type{Val{false}}, ::Type{Val{true}})

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    slave_elements = get_slave_elements(problem)

    function calculate_interface(x::Vector)

        ndofs = round(Int, length(x)/2)
        nnodes = round(Int, ndofs/field_dim)
        u = reshape(x[1:ndofs], field_dim, nnodes)
        la = reshape(x[ndofs+1:end], field_dim, nnodes)
        fc = zeros(u)
        gap = zeros(u)
        C = zeros(la)
        S = Set{Int64}()

        # 1. update nodal normals for slave elements
        Q = [0.0 -1.0; 1.0 0.0]
        normals = zeros(u)
        for element in slave_elements
            conn = get_connectivity(element)
            push!(S, conn...)
            gdofs = get_gdofs(element, field_dim)
            X_el = element("geometry", time)
            u_el = Field(Vector[u[:,i] for i in conn])
            x_el = X_el + u_el
            for ip in get_integration_points(element, 3)
                dN = get_dbasis(element, ip, time)
                N = element(ip, time)
                t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
                normals[:, conn] += ip.weight*Q*t'*N
            end
        end
        for i in 1:size(normals,2)
            normals[:,i] /= norm(normals[:,i])
        end
        # swap element normals in 2d if they point to inside of body
        if props.rotate_normals
            for i=1:size(normals,2)
                normals[:,i] = -normals[:,i]
            end
        end
        normals2 = Dict()
        for j in S
            normals2[j] = normals[:,j]
        end
        update!(slave_elements, "normal", time => normals2)

        # 2. loop all slave elements
        for slave_element in slave_elements

            slave_element_nodes = get_connectivity(slave_element)
            X1 = slave_element("geometry", time)
            u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
            x1 = X1 + u1
            la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
            n1 = Field(Vector[normals[:,i] for i in slave_element_nodes])
            nnodes = size(slave_element, 2)

            # 3. loop all master elements
            for master_element in slave_element("master elements", time)

                master_element_nodes = get_connectivity(master_element)
                X2 = master_element("geometry", time)
                u2 = Field(Vector[u[:,i] for i in master_element_nodes])
                x2 = X2 + u2

                #x1_midpoint = 1/2*(x1[1]+x1[2])
                #x2_midpoint = 1/2*(x2[1]+x2[2])
                #distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint))
                #distance > props.maximum_distance && continue

                # calculate segmentation: we care only about endpoints
                xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
                xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2])
                xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
                l = 1/2*abs(xi1[2]-xi1[1])
                isapprox(l, 0.0) && continue # no contribution in this master element

                De = zeros(nnodes, nnodes)
                Me = zeros(nnodes, nnodes)
                for ip in get_integration_points(slave_element, 3)
                    # jacobian of slave element in deformed state
                    dN = get_dbasis(slave_element, ip, time)
                    j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
                    w = ip.weight*norm(j)*l
                    xi = ip.coords[1]
                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                    N1 = get_basis(slave_element, xi_s, time)
                    De += w*diagm(vec(N1))
                    Me += w*N1'*N1
                end
                Ae = De*inv(Me)

                slave_dofs = get_gdofs(slave_element, field_dim)
                master_dofs = get_gdofs(master_element, field_dim)

                # 4. loop integration points of segment
                for ip in get_integration_points(slave_element, 3)
                    # jacobian of slave element in deformed state
                    dN = get_dbasis(slave_element, ip, time)
                    j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
                    w = ip.weight*norm(j)*l

                    # project gauss point from slave element to master element
                    xi = ip.coords[1]
                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                    N1 = vec(get_basis(slave_element, xi_s, time))
                    x_s = N1*x1 # coordinate in gauss point
                    n_s = N1*n1 # normal direction in gauss point
                    t_s = Q'*n_s # tangent direction in gauss point
                    xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
                    N2 = vec(get_basis(master_element, xi_m, time))
                    x_m = N2*x2
                    Phi = Ae*N1

                    la_s = Phi*la1 # traction force in gauss point
                    gn = -dot(n_s, x_s - x_m) # normal gap

                    fc[:,slave_element_nodes] += w*la_s*N1'
                    fc[:,master_element_nodes] -= w*la_s*N2'
                    gap[1,slave_element_nodes] += w*gn*Phi'
                    #gap[1,slave_element_nodes] += w*gn*N1'

                end # done integrating segment

            end # master elements done

        end # slave elements done

        # at this point we have calculated contact force fc and gap for all slave elements.
        # next task is to find out are they in contact or not and remove inactive nodes

        #nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap))))
        #info("gap: $nzgap")

        for (i, j) in enumerate(sort(collect(S)))
#           if j in props.always_inactive
#               info("special node $j always inactive")
#               C[:,j] = la[:,j]
#               continue
#           end
            n = normals[:,j]
            t = Q'*n
            lan = dot(n, la[:,j])
            lat = dot(t, la[:,j])

            if lan - gap[1, j] > 0
#               info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t))")
                C[1,j] += gap[1, j]
                C[2,j] += lat
            else
                C[:,j] = la[:,j]
            end
        end

        return vec([fc C])

    end

    # x doesn't mean deformed configuration here
    x = [problem.assembly.u; problem.assembly.la]
    if length(x) == 0
        error("2d autodiff contact problem: initialize problem.assembly.u & la before solution")
    end
    A = ForwardDiff.jacobian(calculate_interface, x)
    b = calculate_interface(x)
    A = sparse(A)
    b = sparse(b)
    SparseMatrix.droptol!(A, 1.0e-12)
    SparseMatrix.droptol!(b, 1.0e-12)

    ndofs = round(Int, length(x)/2)
    K = A[1:ndofs,1:ndofs]
    C1 = transpose(A[1:ndofs,ndofs+1:end])
    C2 = A[ndofs+1:end,1:ndofs]
    D = A[ndofs+1:end,ndofs+1:end]
    f = -b[1:ndofs]
    g = -b[ndofs+1:end]

    empty!(problem.assembly)
    add!(problem.assembly.K, K)
    add!(problem.assembly.C1, C1)
    add!(problem.assembly.C2, C2)
    add!(problem.assembly.D, D)
    add!(problem.assembly.f, f)
    add!(problem.assembly.g, g)

end


text:
"""
Banzhaf (1965) value
"""
function fᵦⁱ(G::Game, i::Int)
    N = length(G.N)

    banzhaf = ∑(
        G.v(S) - G.v(setdiff(S, i))
        for S in subsets(G.N)
            if i ∈ S
    )

    return banzhaf ./ (2^N - 1)
end


text:
function maximize(a::NamedArray)

  amat = copy(a)
  #local sr
  #local one, two, onearrow, onearc, twoarrow, twoarc
  na = size(amat, 1)
  at = findall(amat .+ transpose(amat) .+ I(na) .== 0)
  if length(at) > 0
    for i in 1:length(at)
      s = [at[i][1], at[i][2]]
      st = Int[]
      while true
        st = s
        for j in s
          for k in 1:na
            if mod10(amat[k, j]) == 1
              s = vcat(k, filter(x -> !(k in s), s))
            end
          end
        end
        s == st && break
      end
      one = at[i][1]
      two = at[i][2]
      onearrow = Int[]
      onearc = Int[]
      twoarrow = Int[]
      twoarc = Int[]
      sr = s
      #println("sr = $sr")
      #sr = filter(x -> !(at[i][1] in sr), sr)
      #sr = filter(x -> !(at[i][2] in sr), sr)
      #println("sr = $sr")
      if length(sr) > 0
        for j in sr
          if mod10(amat[one, j]) == 1
            onearrow = vcat(onearrow, j)
          end
          if amat[one, j] > 99
            onearc = vcat(onearc, j)
          end
        end
        for j in sr
          if mod10(amat[two, j]) == 1
            for k in onearc
              if j == k || amat[j, k] > 99
                amat[at[i][2], at[i][1]] = 1
              end
              twoarrow = vcat(twoarrow, j)
            end
          end
          if amat[two, j] > 99
            for k in onearrow
              if j == k || amat[j, k] > 99
                amat[at[i][1], at[i][2]] = 1
              end
            end
            for k in onearc
              if j == k || amat[j, k] > 99
                amat[at[i][1], at[i][2]] = 100
                amat[at[i][2], at[i][1]] = 100
              end
            end
            twoarc = vcat(twoarc, j)
          end
        end
      end
      if length(vcat(onearc, onearrow, twoarc, twoarrow)) > 0

        for j in vcat(onearc,onearrow,twoarc,twoarrow)
          sr = filter(x -> !(j in sr), sr)
        end
        onearcn = Int[]
        twoarcn = Int[]
        onearrown = Int[]
        twoarrown = Int[]
        while true
          for l in onearc
            for j in sr
               if amat[l, j] > 99
                onearcn = vcat(onearcn, j)
              end
            end
          end
          for l in onearrow
            for j in sr
              if amat[l, j] > 99
                onearrown = vcat(onearrown, j)
              end
            end
          end
          for l in twoarc
            for j in sr
              if amat[l, j] > 99
                for k in onearrow
                  if j==k || amat[j, k] > 99
                    amat[at[i][1], at[i][2]] = 1
                  end
                end
                for k in onearc
                  if j==k || amat[j, k] > 99
                    amat[at[i][1], at[i][2]] = 100
                    amat[at[i][2], at[i][1]] = 100
                  end
                end
                twoarcn = vcat(twoarcn, j)
              end
            end
          end
          for l in twoarrow
            for j in sr
              if amat[l,j] > 99
                for k in onearc
                  if j==k || amat[j, k] > 99
                    amat[at[i][1], at[i][2]] = 100
                    amat[at[i][2], at[i][1]] = 100
                  end
                end
                twoarrown = vcat(twoarrown, j)
              end
            end
          end

          length(vcat(onearcn,onearrown,twoarcn,twoarrown)) == 0 && break

          for j in vcat(onearcn,onearrown,twoarcn,twoarrown)
            sr = filter(x -> !(j in sr), sr)
          end

          onearc = onearcn
          twoarc = twoarcn
          onearrow = onearrown
          twoarrow = twoarrown
          onearcn = Int[]
          twoarcn = Int[]
          onearrown = Int[]
          twoarrown = Int[]
        end
      end
    end
  end

  amat

end

export
  maximize


text:
struct NNPACKError <: Exception
    code::nnp_status
    msg::AbstractString
end

Base.show(io::IO, err::NNPACKError) = print(io, "NNPACKError(code $(err.code), $(err.msg))")

function NNPACKError(status::nnp_status)
    msg = "NNPACK STATUS SUCCESS"
    if status == nnp_status_invalid_batch_size
        msg = "NNPACK STATUS INVALID BATCH SIZE"
    elseif status == nnp_status_invalid_channels
        msg = "NNPACK STATUS INVALID CHANNELS"
    elseif status == nnp_status_invalid_input_channels
        msg = "NNPACK STATUS INVALID INPUT CHANNELS"
    elseif status == nnp_status_invalid_output_channels
        msg = "NNPACK STATUS INVALID OUTPUT CHANNELS"
    elseif status == nnp_status_invalid_input_size
	    msg = "NNPACK STATUS INVALID INPUT SIZE"
    elseif status == nnp_status_invalid_input_stride
        msg = "NNPACK STATUS INVALID INPUT STRIDE"
    elseif status == nnp_status_invalid_input_padding
        msg = "NNPACK STATUS INVALID INPUT PADDING"
    elseif status == nnp_status_invalid_kernel_size
        msg = "NNPACK STATUS INVALID KERNEL SIZE"
    elseif status == nnp_status_invalid_pooling_size
        msg = "NNPACK STATUS INVALID POOLING SIZE"
    elseif status == nnp_status_invalid_pooling_stride
        msg = "NNPACK STATUS INVALID POOLING STRIDE"
    elseif status == nnp_status_invalid_algorithm
        msg = "NNPACK STATUS INVALID ALGORITHM"
    elseif status == nnp_status_invalid_transform_strategy
        msg = "NNPACK STATUS INVALID TRANSFORM STRATEGY"
    elseif status == nnp_status_invalid_output_subsampling
        msg = "NNPACK STATUS INVALID OUTPUT SUBSAMPLING"
    elseif status == nnp_status_invalid_activation
        msg = "NNPACK STATUS INVALID ACTIVATION"
    elseif status == nnp_status_invalid_activation_parameters
        msg = "NNPACK STATUS INVALID ACTIVATION PARAMETERS"
    elseif status == nnp_status_unsupported_input_size
        msg = "NNPACK STATUS UNSUPPORTED INPUT SIZE"
    elseif status == nnp_status_unsupported_input_stride
        msg = "NNPACK STATUS UNSUPPORTED INPUT STRIDE"
    elseif status == nnp_status_unsupported_input_padding
        msg = "NNPACK STATUS UNSUPPORTED INPUT PADDING"
    elseif status == nnp_status_unsupported_kernel_size
        msg = "NNPACK STATUS UNSUPPORTED KERNEL SIZE"
    elseif status == nnp_status_unsupported_pooling_size
        msg = "NNPACK STATUS UNSUPPORTED POOLING SIZE"
    elseif status == nnp_status_unsupported_pooling_stride
        msg = "NNPACK STATUS UNSUPPORTED POOLING STRIDE"
    elseif status == nnp_status_unsupported_algorithm
        msg = "NNPACK STATUS UNSUPPORTED ALGORITHM"
    elseif status == nnp_status_unsupported_transform_strategy
        msg = "NNPACK STATUS UNSUPPORTED TRANSFORM STRATEGY"
    elseif status == nnp_status_unsupported_activation
        msg = "NNPACK STATUS UNSUPPORTED ACTIVATION"
    elseif status == nnp_status_unsupported_activation_parameters
        msg = "NNPACK STATUS UNSUPPORTED ACTIVATION PARAMETERS"
    elseif status == nnp_status_uninitialized
        msg = "NNPACK STATUS UNINITIALIZED"
    elseif status == nnp_status_unsupported_hardware
        msg = "NNPACK STATUS UNSUPPORTED HARDWARE"
    elseif status == nnp_status_out_of_memory
        msg = "NNPACK STATUS OUT OF MEMORY"
    elseif status == nnp_status_insufficient_buffer
        msg = "NNPACK STATUS INSUFFICIENT BUFFER"
    elseif status == nnp_status_misaligned_buffer
        msg = "NNPACK STATUS MISALIGNED BUFFER"
    end
    NNPACKError(status, msg)
end

macro check(nnp_func)
    quote
        local err::nnp_status
        err = $(esc(nnp_func))
        if err != nnp_status_success
            throw(NNPACKError(err))
        end
        err
    end
end


text:
# issue 505
@testset "mdd+indent" begin
    h = """
        @def x = Dict(
            :a => 5,
            :b => 7
            )
        A
        ```julia:ex
        locvar(:x)[:a]
        ```
        \\show{ex}
        """ |> fd2html
    @test occursin(
        "5", h)
    h = """
        @def x = Dict(
            :a => (1, 2, 3),
            :b => 7
            )
        A
        ```julia:ex
        locvar(:x)[:a][1]
        ```
        \\show{ex}
        """ |> fd2html
    @test occursin(
        "1", h)
    h = """
        @def x = Dict(
            :a => (1,
                   2,
                   3),
            :b => 7
            )
        A
        ```julia:ex
        locvar(:x)[:a][2]
        ```
        \\show{ex}
        """ |> fd2html
    @test occursin(
        "2", h)
end


text:
#
# Edit `Port for Communication with the Julia Process` to 9999 on Atom Juno
#

port = 9999


text:
abstract type BoundingBox end

struct OBB <: BoundingBox
    c::SVector{3,Float64}
    e::SVector{3,Float64}
    R::SMatrix{3,3,Float64,9}
    OBB(c::SVector{3,Float64}, e::SVector{3,Float64}, R::SMatrix{3,3,Float64,9}) = new(c, e, R)
end

function (::Type{BB_Type})(a::BB_Type, b::BB_Type) where {BB_Type <: BoundingBox}
    min_1, max_1 = calc_min_max(a)
    min_2, max_2 = calc_min_max(b)
    return calc_obb(SVector{4,SVector{3,Float64}}(min_1, max_1, min_2, max_2))
end


text:
#xor-ml.jl - set up to test the xor function in julia.

using SigmoidNumbers
using GenML

import GenML.FCL: FullyConnectedTransition, matrixfma, dxasychainrule, scaledsubtract, scaledouterproductfma, reversematrixfma
import GenML.MLP: MultilayerPerceptron
import GenML.dxasy
import GenML.VoidableDeltas
import GenML.BackpropStorage
import GenML.backpropagate!
import GenML.evaluate!

import SigmoidNumbers.Sigmoid
import SigmoidNumbers.__round

GenML.TF.sigmoid(x::Posit) = SigmoidNumbers.pseudologistic(x)
GenML.TF.dasysigmoid(x::Posit) = SigmoidNumbers.delta_psl(x)
GenML.CF.logloss(x::Posit) = SigmoidNumbers.pseudologcost(x)

const clamplist = Symbol[]

#populate the clamplist using values from the command line
__print = false
all = false
count = 100
for arg in ARGS
  if arg == "print"
    __print = true
    count = 1
  end
  if arg == "all"
    all = true
  end
  push!(clamplist, Symbol(arg))
end

#create a "clamping" macro that will clamp values of an array down to a lower posit,
#givn an identifier saying we should clamp this.
macro clamp(value, identifier)

  extracted_symbol = identifier.args[1]

  pstmt = __print ? :(println($value)) : :()

  if all || (extracted_symbol in clamplist)
    esc(quote
      $pstmt
      if isa($value, AbstractArray)
        for idx = 1:length($value)
          if isa($value[idx], Posit)
            temp = reinterpret(Posit{8}, $value[idx])
            temp = SigmoidNumbers.__round(temp)
            $value[idx] = reinterpret(F, temp)
          end
        end
      else
        temp = reinterpret(Posit{8}, $value)
        temp = SigmoidNumbers.__round(temp)
        $value = reinterpret(F, temp)
      end
      $pstmt
    end)
  end
end

__print && println(clamplist)
################################################################################

function GenML.FCL.matrixfma{F<:Posit}(output, mtx, input::AbstractVector{F}, bias, isize, osize)
  for idx = 1:osize
    @inbounds output[idx] = bias[idx]

    @clamp output[idx] :matrixfma_acc

    @inbounds for jdx = 1:isize

      matrixfma_product = mtx[idx, jdx] * input[jdx]

      @clamp matrixfma_product :matrixfma_prod

      output[idx] += matrixfma_product

      @clamp output[idx] :matrixfma_acc
    end
  end
end

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

function GenML.FCL.reversematrixfma{F<:Posit}(input_array::AbstractVector{F}, matrix::AbstractMatrix{F}, output_array::AbstractVector{F}, i, o)
  for idx = 1:i
      input_array[idx] = zero(F)

      @clamp input_array[idx] :reversematrixfma_acc

    for jdx = 1:o

      reversematrixfma_product = matrix[jdx, idx] * output_array[jdx]

      @clamp reversematrixfma_product :reversematrixfma_prod

      input_array[idx] += reversematrixfma_product

      @clamp input_array[idx] :reversematrixfma_acc

    end
  end
end

function GenML.FCL.scaledouterproductfma{F<:Posit}(matrix::AbstractMatrix{F}, output_deltas::AbstractVector{F}, input_array::AbstractVector, alpha::F, i, o)
  for idx = 1:o
    for jdx = 1:i

      scaledouterproductfma_prod = alpha * (output_deltas[idx] * input_array[jdx])

      @clamp scaledouterproductfma_prod :scaledouterproductfma_prod

      matrix[idx,jdx] -= scaledouterproductfma_prod # + lambda * mlp.transitions[l_idx - 1])

      @clamp matrix[idx,jdx] :scaledouterproductfma_acc
    end
  end
end

function GenML.FCL.scaledsubtract{F<:Posit}(target_vector::AbstractVector{F}, value_vector::AbstractVector{F}, alpha::F, o)
  for idx = 1:o
    target_vector[idx] -= alpha * value_vector[idx]

    @clamp target_vector[idx] :scaledsubtract
  end
end

function GenML.FCL.dxasychainrule{F<:Posit}(outer_differential::AbstractVector{F}, inner_value::AbstractVector{F}, f::Function, o)
  for idx = 1:o
    outer_differential[idx] = outer_differential[idx] .* (dxasy(f))(inner_value[idx])

    @clamp outer_differential[idx] :dxasychainrule
  end
end

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

#MONKEY PATCH MLP BACKPROPAGATION ALGO TO HAVE CLAMPS EVERYWHERE
#implementation of the backpropagation algorithm for multilayer perceptrons.
@generated function backpropagate!{F <: Posit, i, o, tf}(fcl::FullyConnectedTransition{F, i, o, tf},
#==#                                            input::AbstractArray,
#==#                                            output::AbstractArray{F},
#==#                                            output_deltas::AbstractArray{F},
#==#                                            input_deltas::VoidableDeltas{F} = nothing)

  code = quote
    #for now.
    alpha = F(0.1)

    @clamp input         :backprop_input
    @clamp output        :backprop_output
    @clamp output_deltas :backprop_output_deltas

    #overwrite the output delta values with the adjusted values taking the
    dxasychainrule(output_deltas, output, tf, o)

    @clamp output_deltas :backprop_output_transition

    scaledsubtract(fcl.bias, output_deltas, alpha, o)

    scaledouterproductfma(fcl.transition, output_deltas, input, alpha, i, o)
  end

  (input_deltas == Void) && return code

  #handle the backpropagating delta data.
  quote
    $code
    reversematrixfma(input_deltas, fcl.transition, output_deltas, i, o)

    @clamp input_deltas :backprop_input_deltas
  end
end

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

function evaluate!{F <: Posit, i, o, tf}(output::AbstractArray{F}, fcl::FullyConnectedTransition{F, i, o, tf}, input::AbstractArray)

  @clamp input :evaluate_input

  matrixfma(output, fcl.transition, input, fcl.bias, i, o)

  @clamp output :evaluate_transition

  for idx = 1:o
    output[idx] = tf(output[idx])
  end

  @clamp output :evaluate_output
end

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

function xortrain(PType)

  println("working on unreliable xor data set with backpropagation")

  input_matrix = rand(Bool, 10, 500)
  training_results = Array{Bool,2}(1,500)
  training_results[:] = [input_matrix[1, col] $ input_matrix[2, col] $ (rand() < 0.05) for col in 1:size(input_matrix,2)]

  xornet = GenML.MLP.MultilayerPerceptron{PType,(10,2,1)}(randn)

  for rounds = 1:25
  for idx = 1:500
    input_column = input_matrix[:, idx]
    training_column = training_results[:, idx]

    GenML.Optimizers.backpropagationoptimize(xornet, input_column, training_column, GenML.CF.crossentropy)
  end
  end

  #verify that the optimization has resulted in a good data set.
  wrongcount = 0
  for x = 1:50
    input_vector = rand(Bool, 10)
    wrongcount += (xornet(input_vector)[1] > PType(0.5)) != (input_vector[1] $ input_vector[2])
  end
  println("incorrect responses for noisy xor: $wrongcount / 50 => $(2 * (50 - wrongcount)) % correct")

  return wrongcount
end

Base.randn{PType <: Sigmoid}(::Type{PType}, dims::Integer...) = map(PType, randn(dims...))

const trials = count

results = zeros(Int, trials)
PType = Posit{16}

for idx = 1:trials
  println("posit $(PType): round $idx")
  results[idx] = xortrain(PType)
end

println("result: ", sum(results) / trials)

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


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

JLLWrappers.@generate_main_file_header("libcleri")
JLLWrappers.@generate_main_file("libcleri", UUID("cdc7adba-bef8-5cba-a7ee-c792dee3081e"))
end  # module libcleri_jll


text:
function _weight_array_to_dict(w::Array{Weight}=Weight[])
    w_dict = Dict()
    for w_dummy::Weight in w
        w_dict[w_dummy.node] = w_dummy.weight
    end
    return w_dict
end

function _label_map_array_to_dict(l::Array{NodeLabel}=NodeLabel[])
    l_dict = Dict()
    l_s_dict = Dict()
    l_i_dict = Dict()
    for l_dummy::NodeLabel in l
        l_dict[l_dummy.node] = l_dummy.label
        l_s_dict[l_dummy.node] = l_dummy.sequence
        l_i_dict[l_dummy.node] = l_dummy.info
    end
    return l_dict, l_s_dict, l_i_dict
end

function _weight_from_dict(w_dict::Dict)
    w::Array{Weight} = Weight[]
    for (k, v) in w_dict
        push!(w, Weight(k, v))
    end
    return w
end

function _label_from_dict(l_dict::Dict, l_s_dict::Dict, l_i_dict::Dict)
    l::Array{NodeLabel} = NodeLabel[]
    for (k, v) in l_dict
        s = l_s_dict[v]
        i = l_i_dict[v]
        push!(l, NodeLabel(k, v, s, i))
    end
    return l
end

function _edge_label_from_dict(e_o_dict::Dict, e_i_dict::Dict)
    e::Array{EdgeLabel} = EdgeLabel[]
    for (k, v) in e_o_dict
        e_dummy = split(k, "::")
        s = e_dummy[1]
        d = e_dummy[2]
        o = v
        i = e_i_dict[k]
        push!(e, EdgeLabel(s, d, o, i))
    end
    return e
end

function _edge_label_array_to_dict(e::Array{EdgeLabel}=EdgeLabel[])
    e_o_dict = Dict()
    e_i_dict = Dict()
    for e_dummy::EdgeLabel in e
        e_name = e_dummy.src * "::" * e_dummy.dst
        e_o_dict[e_name] = e_dummy.overlap
        e_i_dict[e_name] = e_dummy.info
    end
    return e_o_dict, e_i_dict
end

function _sg_map_array_to_dict(sg::Array{SubGraphMap}=SubGraphMap[])
    sg_dict = Dict()
    for sg_dummy::SubGraphMap in sg
        sg_dict[sg_dummy.node] = sg_dummy.master_node
    end
    return sg_dict
end

function _g_opt_to_path(g_opt::SimpleDiGraph, source_nodes::Array{Int64})
    dis_arr = []
    edges_arr = []
    for node in source_nodes
        edges = a_star(g_opt, node, nv(g_opt))
        push!(edges_arr, edges)
        push!(dis_arr, length(edges))
    end
    max_edges = edges_arr[findmax(dis_arr)[2]]
    path = Int64[]
    for e in max_edges
        push!(path, e.src)
    end
    return path
end

function _label_path(graph_result::GraphResult, path::Array{Int64})
    l_dict, l_s_dict, l_i_dict = _label_map_array_to_dict(graph_result.node_label)
    label_path = String[]
    for node in path
        push!(label_path, l_dict[node])
    end
    return label_path
end


text:
using Firebase
using Test
using HTTP, JSON

@testset "Firebase" begin
    include("realtime.jl")
    include("firestore.jl")

end


text:
function nlp_mi_007_020(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - infeasible model, due to integrality

    model = Model(optimizer)

    @variable(model, -2 <= x <= 3, Int)
    @variable(model, y, Bin)

    @NLconstraint(model, (x - 0.5)^2 + (4 * y - 2)^2 <= 3)

    optimize!(model)

    return check_status(
        model,
        INFEASIBLE_PROBLEM,
        termination_target,
        primal_target,
    )
end


text:
### BasicContUnvParameter

# Guidelines for usage of inner constructors of continuous parameter types:
# 1) Function fields have higher priority than implicitly derived definitions via the pdf field
# 2) Target-related fields have higher priority than implicitly derived likelihood+prior fields
# 3) Upto-related fields have higher priority than implicitly derived Function tuples

### BasicContUnvParameter

mutable struct BasicContUnvParameter <: Parameter{Continuous, Univariate}
  key::Symbol
  index::Integer
  pdf::Union{ContinuousUnivariateDistribution, Void}
  prior::Union{ContinuousUnivariateDistribution, Void}
  setpdf::Union{Function, Void}
  setprior::Union{Function, Void}
  loglikelihood!::Union{Function, Void}
  logprior!::Union{Function, Void}
  logtarget!::Union{Function, Void}
  gradloglikelihood!::Union{Function, Void}
  gradlogprior!::Union{Function, Void}
  gradlogtarget!::Union{Function, Void}
  tensorloglikelihood!::Union{Function, Void}
  tensorlogprior!::Union{Function, Void}
  tensorlogtarget!::Union{Function, Void}
  dtensorloglikelihood!::Union{Function, Void}
  dtensorlogprior!::Union{Function, Void}
  dtensorlogtarget!::Union{Function, Void}
  uptogradlogtarget!::Union{Function, Void}
  uptotensorlogtarget!::Union{Function, Void}
  uptodtensorlogtarget!::Union{Function, Void}
  diffmethods::Union{DiffMethods, Void}
  diffopts::Union{DiffOptions, Void}
  states::VariableStateVector
  state::ParameterState{Continuous, Univariate}

  function BasicContUnvParameter(
    key::Symbol,
    index::Integer,
    pdf::Union{ContinuousUnivariateDistribution, Void},
    prior::Union{ContinuousUnivariateDistribution, Void},
    setpdf::Union{Function, Void},
    setprior::Union{Function, Void},
    ll::Union{Function, Void},
    lp::Union{Function, Void},
    lt::Union{Function, Void},
    gll::Union{Function, Void},
    glp::Union{Function, Void},
    glt::Union{Function, Void},
    tll::Union{Function, Void},
    tlp::Union{Function, Void},
    tlt::Union{Function, Void},
    dtll::Union{Function, Void},
    dtlp::Union{Function, Void},
    dtlt::Union{Function, Void},
    uptoglt::Union{Function, Void},
    uptotlt::Union{Function, Void},
    uptodtlt::Union{Function, Void},
    diffmethods::Union{DiffMethods, Void},
    diffopts::Union{DiffOptions, Void},
    states::VariableStateVector,
    state::ParameterState{Continuous, Univariate}
  )
    args = (setpdf, setprior, ll, lp, lt, gll, glp, glt, tll, tlp, tlt, dtll, dtlp, dtlt, uptoglt, uptotlt, uptodtlt)
    fnames = fieldnames(BasicContUnvParameter)[5:21]

    # Check that all generic functions have correct signature
    for i in 1:17
      if isa(args[i], Function) &&
        isa(args[i], Function) &&
        !method_exists(args[i], (BasicContUnvParameterState, VariableStateVector))
        error("$(fnames[i]) has wrong signature")
      end
    end

    new(
      key,
      index,
      pdf,
      prior,
      setpdf,
      setprior,
      ll,
      lp,
      lt,
      gll,
      glp,
      glt,
      tll,
      tlp,
      tlt,
      dtll,
      dtlp,
      dtlt,
      uptoglt,
      uptotlt,
      uptodtlt,
      diffmethods,
      diffopts,
      states,
      state
    )
  end
end

function BasicContUnvParameter!(
  parameter::BasicContUnvParameter,
  setpdf::Union{Function, Void},
  setprior::Union{Function, Void},
  ll::Union{Function, Void},
  lp::Union{Function, Void},
  lt::Union{Function, Void},
  gll::Union{Function, Void},
  glp::Union{Function, Void},
  glt::Union{Function, Void},
  tll::Union{Function, Void},
  tlp::Union{Function, Void},
  tlt::Union{Function, Void},
  dtll::Union{Function, Void},
  dtlp::Union{Function, Void},
  dtlt::Union{Function, Void},
  uptoglt::Union{Function, Void},
  uptotlt::Union{Function, Void},
  uptodtlt::Union{Function, Void}
)
  args = (setpdf, setprior, ll, lp, lt, gll, glp, glt, tll, tlp, tlt, dtll, dtlp, dtlt, uptoglt, uptotlt, uptodtlt)

  # Define setpdf (i = 1) and setprior (i = 2)
  for (i, setter, distribution) in ((1, :setpdf, :pdf), (2, :setprior, :prior))
    setfield!(
      parameter,
      setter,
      if isa(args[i], Function)
        (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          setfield!(parameter, distribution, args[i](_state, _states))
      else
        nothing
      end
    )
  end

  # Define loglikelihood! (i = 3) and gradloglikelihood! (i = 6)
  # plfield stands for parameter likelihood-related field respectively
  for (i, plfield) in ((3, :loglikelihood!), (6, :gradloglikelihood!))
    setfield!(
      parameter,
      plfield,
      if isa(args[i], Function)
        _state::BasicContUnvParameterState -> args[i](_state, parameter.states)
      else
        nothing
      end
    )
  end

  # Define logprior! (i = 4) and gradlogprior! (i = 7)
  # ppfield and spfield stand for parameter prior-related field and state prior-related field repsectively
  for (i , ppfield, spfield, f) in ((4, :logprior!, :logprior, logpdf), (7, :gradlogprior!, :gradlogprior, gradlogpdf))
    setfield!(
    parameter,
      ppfield,
      if isa(args[i], Function)
        _state::BasicContUnvParameterState -> args[i](_state, parameter.states)
      else
        if (
            isa(parameter.prior, ContinuousUnivariateDistribution) &&
            method_exists(f, (typeof(parameter.prior), eltype(parameter.prior)))
          ) ||
          isa(args[2], Function)
          _state::BasicContUnvParameterState -> setfield!(_state, spfield, f(getfield(parameter, :prior), _state.value))
        else
          nothing
        end
      end
    )
  end

  # Define logtarget! (i = 5) and gradlogtarget! (i = 8)
  # ptfield, plfield and ppfield stand for parameter target, likelihood and prior-related field respectively
  # stfield, slfield and spfield stand for state target, likelihood and prior-related field respectively
  for (i , ptfield, plfield, ppfield, stfield, slfield, spfield, f) in (
    (5, :logtarget!, :loglikelihood!, :logprior!, :logtarget, :loglikelihood, :logprior, logpdf),
    (8, :gradlogtarget!, :gradloglikelihood!, :gradlogprior!, :gradlogtarget, :gradloglikelihood, :gradlogprior, gradlogpdf)
  )
    setfield!(
      parameter,
      ptfield,
      if isa(args[i], Function)
        _state::BasicContUnvParameterState -> args[i](_state, parameter.states)
      else
        if isa(args[i-2], Function) && isa(getfield(parameter, ppfield), Function)
          function (_state::BasicContUnvParameterState)
            getfield(parameter, plfield)(_state)
            getfield(parameter, ppfield)(_state)
            setfield!(_state, stfield, getfield(_state, slfield)+getfield(_state, spfield))
          end
        elseif (
            isa(parameter.pdf, ContinuousUnivariateDistribution) &&
            method_exists(f, (typeof(parameter.pdf), eltype(parameter.pdf)))
          ) ||
          isa(args[1], Function)
          _state::BasicContUnvParameterState -> setfield!(_state, stfield, f(getfield(parameter, :pdf), _state.value))
        else
          nothing
        end
      end
    )
  end

  # Define tensorloglikelihood! (i = 9) and dtensorloglikelihood! (i = 12)
  # plfield stands for parameter likelihood-related field respectively
  for (i, plfield) in ((9, :tensorloglikelihood!), (12, :dtensorloglikelihood!))
    setfield!(
      parameter,
      plfield,
      if isa(args[i], Function)
        _state::BasicContUnvParameterState -> args[i](_state, parameter.states)
      else
        nothing
      end
    )
  end

  # Define tensorlogprior! (i = 10) and dtensorlogprior! (i = 13)
  # ppfield stands for parameter prior-related field respectively
  for (i, ppfield) in ((10, :tensorlogprior!), (13, :dtensorlogprior!))
    setfield!(
      parameter,
      ppfield,
      if isa(args[i], Function)
        _state::BasicContUnvParameterState -> args[i](_state, parameter.states)
      else
        nothing
      end
    )
  end

  # Define tensorlogtarget! (i = 11) and dtensorlogtarget! (i = 14)
  for (i , ptfield, plfield, ppfield, stfield, slfield, spfield) in (
    (
      11,
      :tensorlogtarget!, :tensorloglikelihood!, :tensorlogprior!,
      :tensorlogtarget, :tensorloglikelihood, :tensorlogprior
    ),
    (
      14,
      :dtensorlogtarget!, :dtensorloglikelihood!, :dtensorlogprior!,
      :dtensorlogtarget, :dtensorloglikelihood, :dtensorlogprior
    )
  )
    setfield!(
      parameter,
      ptfield,
      if isa(args[i], Function)
        _state::BasicContUnvParameterState -> args[i](_state, parameter.states)
      else
        if isa(args[i-2], Function) && isa(args[i-1], Function)
          function (_state::BasicContUnvParameterState)
            getfield(parameter, plfield)(_state)
            getfield(parameter, ppfield)(_state)
            setfield!(_state, stfield, getfield(_state, slfield)+getfield(_state, spfield))
          end
        else
          nothing
        end
      end
    )
  end

  # Define uptogradlogtarget!
  setfield!(
    parameter,
    :uptogradlogtarget!,
    if isa(args[15], Function)
      _state::BasicContUnvParameterState -> args[15](_state, parameter.states)
    else
      if isa(parameter.logtarget!, Function) && isa(parameter.gradlogtarget!, Function)
        function (state::BasicContUnvParameterState)
          parameter.logtarget!(state)
          parameter.gradlogtarget!(state)
        end
      else
        nothing
      end
    end
  )

  # Define uptotensorlogtarget!
  setfield!(
    parameter,
    :uptotensorlogtarget!,
    if isa(args[16], Function)
      _state::BasicContUnvParameterState -> args[16](_state, parameter.states)
    else
      if isa(parameter.logtarget!, Function) &&
        isa(parameter.gradlogtarget!, Function) &&
        isa(parameter.tensorlogtarget!, Function)
        function (state::BasicContUnvParameterState)
          parameter.logtarget!(state)
          parameter.gradlogtarget!(state)
          parameter.tensorlogtarget!(state)
        end
      else
        nothing
      end
    end
  )

  # Define uptodtensorlogtarget!
  setfield!(
    parameter,
    :uptodtensorlogtarget!,
    if isa(args[17], Function)
      _state::BasicContUnvParameterState -> args[17](_state, parameter.states)
    else
      if isa(parameter.logtarget!, Function) &&
        isa(parameter.gradlogtarget!, Function) &&
        isa(parameter.tensorlogtarget!, Function) &&
        isa(parameter.dtensorlogtarget!, Function)
        function (state::BasicContUnvParameterState)
          parameter.logtarget!(state)
          parameter.gradlogtarget!(state)
          parameter.tensorlogtarget!(state)
          parameter.dtensorlogtarget!(state)
        end
      else
        nothing
      end
    end
  )
end

BasicContUnvParameter(key::Symbol, index::Integer=0; signature::Symbol=:high, args...) =
  BasicContUnvParameter(key, Val{signature}, index; args...)

function BasicContUnvParameter(
  key::Symbol,
  ::Type{Val{:low}},
  index::Integer=0;
  pdf::Union{ContinuousUnivariateDistribution, Void}=nothing,
  prior::Union{ContinuousUnivariateDistribution, Void}=nothing,
  setpdf::Union{Function, Void}=nothing,
  setprior::Union{Function, Void}=nothing,
  loglikelihood::Union{Function, Void}=nothing,
  logprior::Union{Function, Void}=nothing,
  logtarget::Union{Function, Void}=nothing,
  gradloglikelihood::Union{Function, Void}=nothing,
  gradlogprior::Union{Function, Void}=nothing,
  gradlogtarget::Union{Function, Void}=nothing,
  tensorloglikelihood::Union{Function, Void}=nothing,
  tensorlogprior::Union{Function, Void}=nothing,
  tensorlogtarget::Union{Function, Void}=nothing,
  dtensorloglikelihood::Union{Function, Void}=nothing,
  dtensorlogprior::Union{Function, Void}=nothing,
  dtensorlogtarget::Union{Function, Void}=nothing,
  uptogradlogtarget::Union{Function, Void}=nothing,
  uptotensorlogtarget::Union{Function, Void}=nothing,
  uptodtensorlogtarget::Union{Function, Void}=nothing,
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  states::VariableStateVector=VariableState[],
  state::ParameterState{Continuous, Univariate}=BasicContUnvParameterState()
)
  parameter = BasicContUnvParameter(key, index, pdf, prior, fill(nothing, 17)..., diffmethods, diffopts, states, state)

  BasicContUnvParameter!(
    parameter,
    setpdf,
    setprior,
    loglikelihood,
    logprior,
    logtarget,
    gradloglikelihood,
    gradlogprior,
    gradlogtarget,
    tensorloglikelihood,
    tensorlogprior,
    tensorlogtarget,
    dtensorloglikelihood,
    dtensorlogprior,
    dtensorlogtarget,
    uptogradlogtarget,
    uptotensorlogtarget,
    uptodtensorlogtarget
  )

  parameter
end

function BasicContUnvParameter(
  key::Symbol,
  ::Type{Val{:high}},
  index::Integer=0;
  pdf::Union{ContinuousUnivariateDistribution, Void}=nothing,
  prior::Union{ContinuousUnivariateDistribution, Void}=nothing,
  setpdf::Union{Function, Void}=nothing,
  setprior::Union{Function, Void}=nothing,
  loglikelihood::Union{Function, Expr, Void}=nothing,
  logprior::Union{Function, Expr, Void}=nothing,
  logtarget::Union{Function, Expr, Void}=nothing,
  gradloglikelihood::Union{Function, Void}=nothing,
  gradlogprior::Union{Function, Void}=nothing,
  gradlogtarget::Union{Function, Void}=nothing,
  tensorloglikelihood::Union{Function, Void}=nothing,
  tensorlogprior::Union{Function, Void}=nothing,
  tensorlogtarget::Union{Function, Void}=nothing,
  dtensorloglikelihood::Union{Function, Void}=nothing,
  dtensorlogprior::Union{Function, Void}=nothing,
  dtensorlogtarget::Union{Function, Void}=nothing,
  uptogradlogtarget::Union{Function, Void}=nothing,
  uptotensorlogtarget::Union{Function, Void}=nothing,
  uptodtensorlogtarget::Union{Function, Void}=nothing,
  nkeys::Integer=0,
  vfarg::Bool=false,
  diffopts::Union{DiffOptions, Void}=nothing,
  states::VariableStateVector=VariableState[],
  state::ParameterState{Continuous, Univariate}=BasicContUnvParameterState()
)
  inargs = (
    setpdf,
    setprior,
    loglikelihood,
    logprior,
    logtarget,
    gradloglikelihood,
    gradlogprior,
    gradlogtarget,
    tensorloglikelihood,
    tensorlogprior,
    tensorlogtarget,
    dtensorloglikelihood,
    dtensorlogprior,
    dtensorlogtarget,
    uptogradlogtarget,
    uptotensorlogtarget,
    uptodtensorlogtarget
  )

  fnames = fieldnames(BasicContUnvParameterState)[2:13]

  if vfarg
    if nkeys > 0
      if diffopts != nothing
        error("In the case of autodiff, if nkeys is not 0, then vfarg must be false")
      end
    else
      error("If vfarg=$vfarg, nkeys must be positive, got $nkeys")
    end
  else
    @assert nkeys >= 0 "If vfarg=$vfarg, nkeys must be non-negative, got $nkeys"
  end

  if diffopts != nothing
    if diffopts.mode != :forward
      error("Only :forward mode autodiff is available for univariate parameters, got $autodiff mode")
    end

    if diffopts.order != 1
      error("Derivative order must be 1, got order=$order")
    end
  end

  if diffopts != nothing
    for i in 3:5
      if isa(inargs[i], Function) && !isa(inargs[i+3], Function)
        diffopts.targets[i-2] = true
      end
    end
  end

  parameter = BasicContUnvParameter(key, index, pdf, prior, fill(nothing, 17)..., DiffMethods(), diffopts, states, state)

  outargs = Union{Function, Void}[nothing for i in 1:17]

  if vfarg
    for i in 1:2
      if isa(inargs[i], Function)
        outargs[i] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          inargs[i](Any[s.value for s in _states])
      end
    end

    for i in 3:14
      if isa(inargs[i], Function)
        outargs[i] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          setfield!(_state, fnames[i-2], inargs[i](Any[s.value for s in _states]))
      end
    end

    if isa(inargs[15], Function)
      outargs[15] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
        (_state.logtarget, _state.gradlogtarget) = inargs[15](Any[s.value for s in _states])
    end

    if isa(inargs[16], Function)
      outargs[16] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
        (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget) = inargs[16](Any[s.value for s in _states])
    end

    if isa(inargs[17], Function)
      outargs[17] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
        (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget, _state.dtensorlogtarget) =
          inargs[17](Any[s.value for s in _states])
    end
  else
    if nkeys > 0
      for i in 1:2
        if isa(inargs[i], Function)
          outargs[i] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
            inargs[i](_state.value, Any[s.value for s in _states])
        end
      end

      for i in 3:14
        if isa(inargs[i], Function)
          outargs[i] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
            setfield!(_state, fnames[i-2], inargs[i](_state.value, Any[s.value for s in _states]))
        end
      end

      if isa(inargs[15], Function)
        outargs[15] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget) = inargs[15](_state.value, Any[s.value for s in _states])
      end

      if isa(inargs[16], Function)
        outargs[16] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget) =
            inargs[16](_state.value, Any[s.value for s in _states])
      end

      if isa(inargs[17], Function)
        outargs[17] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget, _state.dtensorlogtarget) =
            inargs[17](_state.value, Any[s.value for s in _states])
      end
    else
      for i in 1:2
        if isa(inargs[i], Function)
          outargs[i] = (_state::BasicContUnvParameterState, _states::VariableStateVector) -> inargs[i](_state.value)
        end
      end

      for i in 3:14
        if isa(inargs[i], Function)
          outargs[i] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
            setfield!(_state, fnames[i-2], inargs[i](_state.value))
        end
      end

      if isa(inargs[15], Function)
        outargs[15] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget) = inargs[15](_state.value)
      end

      if isa(inargs[16], Function)
        outargs[16] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget) = inargs[16](_state.value)
      end

      if isa(inargs[17], Function)
        outargs[17] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          (_state.logtarget, _state.gradlogtarget, _state.tensorlogtarget, _state.dtensorlogtarget) =
            inargs[17](_state.value)
      end
    end
  end

  for (i, j, field, distribution, f) in ((4, 2, :logprior, :prior, logpdf), (5, 1, :logtarget, :pdf, logpdf))
    if !isa(inargs[i], Function) &&
      (
        (
          isa(getfield(parameter, distribution), ContinuousUnivariateDistribution) &&
          method_exists(f, (typeof(getfield(parameter, distribution)), eltype(getfield(parameter, distribution))))
        ) ||
        isa(inargs[j], Function)
      )
      outargs[i] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
        setfield!(_state, field, f(getfield(parameter, distribution), _state.value))
    end
  end

  if diffopts != nothing
    for (i, field) in ((3, :closurell), (4, :closurelp), (5, :closurelt))
      if isa(inargs[i], Function)
        setfield!(
          parameter.diffmethods,
          field,
          nkeys == 0 ? inargs[i] : _x::Real -> inargs[i](_x, Any[s.value for s in parameter.states])
        )
      end
    end

    for (i, returnname, diffresult, diffmethod) in (
      (6, :gradloglikelihood, :resultll, :closurell),
      (7, :gradlogprior, :resultlp, :closurelp),
      (8, :gradlogtarget, :resultlt, :closurelt)
    )
      if !isa(inargs[i], Function) && isa(inargs[i-3], Function)
        outargs[i] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
          setfield!(
            _state,
            returnname,
            forward_autodiff_derivative(
              getfield(_state.diffstate, diffresult), getfield(_state.diffmethods, diffmethod), _state.value
            )
          )
      end
    end

    if !isa(inargs[15], Function) && isa(inargs[5], Function)
      outargs[15] = (_state::BasicContUnvParameterState, _states::VariableStateVector) ->
        (_state.logtarget, _state.gradlogtarget) =
          forward_autodiff_upto_derivative(_state.diffstate.resultlt, _state.diffmethods.closurelt, _state.value)
    end
  end

  BasicContUnvParameter!(parameter, outargs...)

  parameter
end

value_support(::Type{BasicContUnvParameter}) = Continuous
value_support(::BasicContUnvParameter) = Continuous

variate_form(::Type{BasicContUnvParameter}) = Univariate
variate_form(::BasicContUnvParameter) = Univariate

default_state_type(::BasicContUnvParameter) = BasicContUnvParameterState

default_state(variable::BasicContUnvParameter, value::N, outopts::Dict) where {N<:Real} =
  BasicContUnvParameterState(
    value,
    (haskey(outopts, :diagnostics) && in(:accept, outopts[:diagnostics])) ? [:accept] : Symbol[],
    variable.diffmethods,
    variable.diffopts
  )

show(io::IO, ::Type{BasicContUnvParameter}) = print(io, "BasicContUnvParameter")


text:
using ForwardDiff
using LinearAlgebra

function marquardt(f, x₀; ε=1e-6, maxiter=1000, λ=10.0^4, α=2)
    x = x₀
    fx = f(x)
    for i in 1:maxiter
        g = ForwardDiff.gradient(f, x)
        norm(g) ≤ ε && return (x=x, converged=true, iters=i)
        x′ = x .- (ForwardDiff.hessian(f, x) + λ*I) \ g
        fx′ = f(x′)
        if fx′ < fx
            λ *= 0.5
            fx = fx′
            x = x′
        else
            λ *= 2.0
        end
    end
    (x=x, converged=false, iters=maxiter)
end

max_stars_repo_path: Chapter04/08. Writing a simple optimization routine/marquardt.jl


text:
"""
    Parameters1PL <: AbstractParametersBinary

Contains information of a set of item parameters (values, bounds, priors, posteiors, chains, expected Fisher information, calibrated) under the 1-parameter logistic model.
"""
mutable struct Parameters1PL <: AbstractParametersBinary
    b::Float64
    bounds_b::Vector{Float64}
    prior::Distributions.UnivariateDistribution
    posterior::Distributions.UnivariateDistribution
    chain::Vector{Float64}
    expected_information::Float64
    calibrated::Bool

    Parameters1PL(b, bounds_b, prior, posterior, chain, expected_information, calibrated) =
        new(b, bounds_b, prior, posterior, chain, expected_information, calibrated)

    # Random Initializers
    function Parameters1PL(dist::Distributions.UnivariateDistribution, bounds_b)
        pars = truncate_rand(dist, bounds_b)
        new(pars[1][1], bounds_b, dist, dist, [pars[1]], 1.0, false)
    end

    function Parameters1PL()
        Parameters1PL(Distributions.Normal(0, 1), [-6.0, 6.0], false)
    end
end

"""
    _empty_chain!(parameters::Parameters1PL)
"""
function _empty_chain!(parameters::Parameters1PL)
    parameters.chain = Vector{Float64}(undef, 0)
end


text:
#=
test_types:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-09
=#

module TestTypes
    # packages
    using Test

    # external packages
    include("composite_types.jl")
    include("union_types.jl")
    include("parametric_types.jl")

    # methods
    function test()
        @testset "Type" begin
            @testset "ComposTypes" begin
                @test_nowarn ComposTypes.main()
            end
            @testset "UnionTypes" begin
                @test_nowarn UnionTypes.main()
            end
            @testset "Parametric" begin
                @test_nowarn Parametric.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestTypes
    TestTypes.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using MinFEM
using LinearAlgebra

function semilinear(mesh::Mesh, L::AbstractMatrix, M::AbstractMatrix,
                    s::AbstractVector, BoundaryIndices::Set{Int64}=[], tol=1e-10)

  y = zeros(mesh.nnodes)

  pde = PDESystem(A=L, b=M*s, bc=zeros(mesh.nnodes), DI=BoundaryIndices)

  res = Inf
  while res > tol
    pde.A = L + asmCubicDerivativeMatrix(mesh, y)
    pde.b = -L*y + M*s - asmCubicTerm(mesh, y)
    refresh(pde)
    solve(pde)

    y += pde.state
    res = norm(pde.state)
  end
  return y
end

function solve_semilinear(file_name::String)

  mesh = import_mesh(file_name)

  L = asmLaplacian(mesh)
  M = asmMassMatrix(mesh)

  f(x) = 100.0*x[1]*x[2]
  s = evaluateMeshFunction(mesh, f)

  boundary = union!(mesh.Boundaries[1001].Nodes,
                    mesh.Boundaries[1002].Nodes,
                    mesh.Boundaries[1003].Nodes,
                    mesh.Boundaries[1004].Nodes);

  y = semilinear(mesh, L, M, s, boundary);
  return true
end

@test solve_semilinear("test_v4.msh")


text:
@testset "Contrast Ingredients" begin
    @test ContrastIngredient(0) == ContrastIngredient("IODINE")
    @test ContrastIngredient(1) == ContrastIngredient("GADOLINIUM")
    @test ContrastIngredient(2) == ContrastIngredient("CARBON")
    @test ContrastIngredient(3) == ContrastIngredient("DIOXIDE")
    @test ContrastIngredient(4) == ContrastIngredient("BARIUM")
    @test ContrastIngredient(5) == ContrastIngredient("XENON")
    @test ContrastIngredient(6) == ContrastIngredient("UnkownContrast")
    @test Symbol(String(ContrastIngredient(6))) == Symbol(ContrastIngredient("PLUTIONIUM"))
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the LaTeX backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Issue #38
# ==============================================================================

@testset "Issue #38 - isoverlong is not defined" begin
    data = ["2.0 ± 1" "3.0 ± 1"]

    expected = """
\\begin{table}
  \\begin{tabular}{rr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} \\\\\\hline
    2.0 ± 1 & 3.0 ± 1 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data, backend = :latex)
    @test result == expected
end


text:
# Functions that shift the indexes of the points in an inverted index
function _shift_up_inverse_index!(inverse_index::InvertedIndex{I,U}, by::I=one(I)) where{I,U}
    @inbounds @simd for i = eachindex(inverse_index)
        inverse_index[i].idxs .+= by
    end
end

function _shift_down_inverse_index!(inverse_index::InvertedIndex{I,U}, by::I=one(I)) where{I,U}
    @inbounds @simd for i = eachindex(inverse_index)
        inverse_index[i].idxs .-= by
    end
end

function _shift_inverse_index!(inverse_index::InvertedIndex{I,U}, point::I, by::I=one(I)) where{I,U}
    @inbounds @simd for i = eachindex(inverse_index)
        inverse_index[i].idxs[inverse_index[i].idxs .> point] .-= by
    end
end

"""
    pop!(ivfadc)

Pops from the index `ivfadc` the point with the highest index and returns it updating
the index as well.
"""
pop!(ivfadc::IVFADCIndex) = _pop!(ivfadc, :last)

"""
    popfirst!(ivfadc)

Pops from the index `ivfadc` the first point and returns it updating the index as well.
"""
popfirst!(ivfadc::IVFADCIndex) = _pop!(ivfadc, :first)

# Utility function for poping
function _pop!(ivfadc::IVFADCIndex{U,I,Dc,Dr,T,Q},
               position::Symbol=:last
              ) where{U,I,Dc,Dr,T,Q}
    @assert length(ivfadc) > 0 "Cannot pop element from empty index"
    cluster = 0  # cluster with max index
    idx = 0      # index in inverted list
    (vecid, shift) = I.(ifelse(position==:last, (length(ivfadc)-1, 0), (0, 1)))
    local point_codes
    @inbounds for i = eachindex(ivfadc.inverse_index)
        if vecid in ivfadc.inverse_index[i].idxs
            cluster = i
            idx = findfirst(isequal(vecid), ivfadc.inverse_index[i].idxs)
            point_codes = ivfadc.inverse_index[i].codes[idx]
        end
    end

    # Get point
    reconstructed = _get_quantizer_vector(ivfadc.coarse_quantizer, cluster) +
                        _decode_point(ivfadc.residual_quantizer, point_codes)

    # Delete index data
    deleteat!(ivfadc.inverse_index[cluster].idxs, idx)
    deleteat!(ivfadc.inverse_index[cluster].codes, idx)

    # Shift index
    _shift_down_inverse_index!(ivfadc.inverse_index, shift)
    return reconstructed
end

function _decode_point(quantizer::QuantizedArrays.OrthogonalQuantizer{U,Dr,T,2},
                       codes::Vector{U}
                      ) where{U,Dr,T}
    m, n = length(quantizer.codebooks), quantizer.dims[1]
    residual = Vector{T}(undef, n)
    @inbounds for i in 1:m
        rr = QuantizedArrays.rowrange(n, m, i)
        residual[rr] .= quantizer.codebooks[i][codes[i]]
    end
    return residual
end

"""
    delete_from_index!(ivfadc, points)

Deletes the points with indices contained in `points` from
the index `ivfadc`.
"""
function delete_from_index!(ivfadc::IVFADCIndex{U,I,Dc,Dr,T,Q},
                            points::Vector{<:Integer}
                           ) where{U,I,Dc,Dr,T,Q}
    shifted_points = I.(points .- 1)  # shift points
    for point in sort(unique(shifted_points), rev=true)
        @inbounds for i = eachindex(ivfadc.inverse_index)
            if point in ivfadc.inverse_index[i].idxs
                pidx = findfirst(isequal(point), ivfadc.inverse_index[i].idxs)
                deleteat!(ivfadc.inverse_index[i].idxs, pidx)
                deleteat!(ivfadc.inverse_index[i].codes, pidx)
                _shift_inverse_index!(ivfadc.inverse_index, point)
                break
            end
        end
    end
end

"""
    push!(ivfadc, point)

Pushes `point` to the end of index `ivfadc`; the point is assigned to a cluster
and its quantized code added to the inverted list corresponding to the cluster.
"""
push!(ivfadc::IVFADCIndex, point) = _push!(ivfadc, point, :last)

"""
    pushfirst!(ivfadc, point)

Pushes `point` to the beginning of index `ivfadc`; the point is assigned to a cluster
and its quantized code added to the inverted list corresponding to the cluster.
"""
pushfirst!(ivfadc::IVFADCIndex, point) = _push!(ivfadc, point, :first)

# Utility function for pushing
function _push!(ivfadc::IVFADCIndex{U,I,Dc,Dr,T,Q},
                point::Vector{T},
                position::Symbol
               ) where{U,I,Dc,Dr,T,Q}
    # Checks and initializations
    nrows, nvectors = size(ivfadc)
    @assert nrows == length(point) "Adding to index requires $nrows-element vectors"
    @assert QuantizedArrays.TYPE_TO_BITS[I] >=
        log2(nvectors+1) "Cannot index, exceeding index capacity of $(Int(typemax(I)+1)) points"

    qpoint, mincluster = _encode_point(ivfadc, point)

    # Insert in the inverted list corresponding to the cluster
    (vecid, shift) = ifelse(position == :first, (0, one(I)), (nvectors, zero(I)))
    _shift_up_inverse_index!(ivfadc.inverse_index, shift)
    push!(ivfadc.inverse_index[mincluster].idxs, vecid)
    push!(ivfadc.inverse_index[mincluster].codes, qpoint)
    return nothing
end

function _encode_point(ivfadc::IVFADCIndex{U,I,Dc,Dr,T,Q},
                       point::Vector{T}
                      ) where{U,I,Dc,Dr,T,Q}
    nrows, nvectors = size(ivfadc)

    # Find belonging cluster (first [1] for (idxs,dist), second for idxs array)
    mincluster = coarse_search(ivfadc.coarse_quantizer, point, 1)[1][1]

    # Quantize residual
    residual = point - _get_quantizer_vector(ivfadc.coarse_quantizer, mincluster)
    quantized_point = vec(QuantizedArrays.quantize_data(ivfadc.residual_quantizer,
                            reshape(residual, nrows, 1)))
    return quantized_point, mincluster
end


text:
# quiver

#JL Plots
Using Plots
n = 7
qf(x,y) = 1 / hypot(x, y)
x = repeat(-3:(2*3)/n:3, 1, n) |> vec
y = repeat(-3:(2*3)/n:3, 1, n)' |> vec
vx = qf.(x,y) .* cos.(atan.(y,x)) |> vec
vy = qf.(x,y) .* sin.(atan.(y,x)) |> vec

Plots.quiver(x, y, quiver=(vx, vy), aspect_ratio=:equal, c=:steelblue, linewidth=2)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
include("../../../iron_line_profile.jl")

mass_vars = [10, 100, 10e6]
spin_vars = [0.0, 0.5, 0.998]
angle_vars = [15.0, 40.0, 85.0]

mass_fixed = 10
spin_fixed = 0.998
angle_fixed = 40.0

# for mass in mass_vars
plt, hmap, title, cache = iron_line_profile(;
                                            mass=mass_fixed,
                                            spin=spin_fixed,
                                            obs_angle=angle_fixed,
                                            tolerance=1e-12,
                                            dtmax=0.5,
                                            size_multiplier=10,
                                            resolution=3000,
                                            output="plot",
                                            nbins=200,
                                            fov=10
                                            )
title!(title, titlefontsize=10)
png(plt, "iron_line_basic.png")
# end

# for spin in spin_vars
#     hmap, cache, title, new_img = iron_line_profile(;
#                                                     spin=spin,
#                                                     mass=mass_fixed,
#                                                     obs_angle=angle_fixed,
#                                                     tolerance=1e-12,
#                                                     dtmax=0.5,
#                                                     size_multiplier=10,
#                                                     resolution=3000,
#                                                     output="plot",
#                                                     nbins=300,
#                                                     fov=10
#                                                     )
#     title!(title)
#     png(hmap, "iron_line_spin=$spin.png")
# end

# for angle in angle_vars
#     hmap, cache, title, new_img = iron_line_profile(;
#                                                     obs_angle=angle,
#                                                     spin=spin_fixed,
#                                                     mass=mass_fixed,
#                                                     tolerance=1e-12,
#                                                     dtmax=0.5,
#                                                     size_multiplier=10,
#                                                     resolution=3000,
#                                                     output="plot",
#                                                     nbins=300,
#                                                     fov=10
#                                                     )
#     title!(title)
#     png(hmap, "iron_line_angle=$angle.png")
# end


text:
using Test

@testset "$(file)" for file in readdir(@__DIR__)
    if file == "runtests.jl" || !endswith(file, ".jl")
        continue
    end
    include(file)
end


text:
using PackageCompiler

target = joinpath("build", "sys.so")

if ispath(target)
    mv(target, target * ".bak"; force = true)
end

imgfile = compile_package("DifferentialEquations", "HypergeometricFunctions",
                          "Interpolations", "Roots")
filename = splitdir(imgfile)[2]
@assert filename == splitdir(target)[2]
mv(imgfile, target)


text:
# PDE

module pde #jl

using ReachabilityModels, MAT

file = matopen(@relpath "pde.mat")

## system matrix
A = float(read(file, "A"))  # the matrix has Int entries

## input matrix
B = read(file, "B")

## state domain
X = Universe(84)

## input domain
U = BallInf([0.75], .25)

function model(X0)
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)
    return IVP(S, X0)
end

end  #jl

# Reachability settings

# Results

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
@testset "InvertDict" begin
    @test invertDict(OrderedDict()) == OrderedDict()
    @test invertDict(OrderedDict(0 => 1, 2 => 3)) == OrderedDict(1 => [0], 3 => [2])
    @test invertDict(OrderedDict(0 => 1, 2 => 3, 3 => 3)) == OrderedDict(1 => [0], 3 => [2])
end

@testset "Labels" begin
    @test sort(collect(values(labelsMap(cl1)))) == UInt32.(1:length(cl1))
    @test sort(collect(values(labelsMap(cl2)))) == UInt32.([1,1,2,2,3,3,4,5,6])
end


text:
using XFloats
using LinearAlgebra
using BenchmarkTools
const BT=BenchmarkTools.DEFAULT_PARAMETERS;
BT.overhead = BenchmarkTools.estimate_overhead(); BT.evals=1; BT.time_tolerance = 2.0e-9; BT.samples = 300;

include("noelide.jl")

include("arithmetic.jl")


text:
@doc raw"""
`PSDE`: Stratonovich Partitioned Stochastic Differential Equation

Defines a partitioned stochastic differential initial value problem
```math
\begin{align*}
\dq (t) &= v(t, q(t)) \, dt + B(t, q(t)) \circ dW , & q(t_{0}) &= q_{0} ,
\dp (t) &= f(t, q(t)) \, dt + G(t, q(t)) \circ dW , & p(t_{0}) &= p_{0}
\end{align*}
```
with the drift vector fields ``v`` and ``f``, diffusion matrices ``B`` and ``G``,
initial conditions ``q_{0}`` and ``p_{0}``, the dynamical variables ``(q,p)`` taking
values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}``, and the m-dimensional Wiener process W

### Fields

* `d`:  dimension of dynamical variable ``q`` and the vector field ``v``
* `m`:  dimension of the Wiener process
* `ni`: number of initial conditions
* `ns`: number of sample paths
* `v`:  function computing the drift vector field for the position variable q
* `f`:  function computing the drift vector field for the momentum variable p
* `B`:  function computing the d x m diffusion matrix for the position variable q
* `G`:  function computing the d x m diffusion matrix for the momentum variable p
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q`` (may be a random variable itself)
* `p₀`: initial condition for dynamical variable ``p`` (may be a random variable itself)

The functions `v`, `f`, 'B' and `G`, providing the drift vector fields and diffusion matrices, take four arguments,
`v(t, q, p, v)`, `f(t, q, p, f)`, `B(t, q, p,  B)` and `G(t, q, p, G)`, where `t` is the current time, `(q, p)` is the
current solution vector, and `v`, `f`, 'B' and `G` are the variables which hold the result
of evaluating the vector fields ``v``, ``f`` and the matrices ``B``, ``G`` on `t` and `(q,p)`.

### Example

```julia
    function v(λ, t, q, v)
        v[1] = λ*q[1]
        v[2] = λ*q[2]
    end

    function B(μ, t, q, B)
        B[1] = μ*q[1]
        B[2] = μ*q[2]
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ  = 2.
    μ  = 1.

    v_sde = (t, q, v) -> v(λ, t, q, v)
    B_sde = (t, q, B) -> B(μ, t, q, B)

    sde = SDE(v_sde, B_sde, t₀, q₀)
```
"""
struct PSDE{dType <: Number, tType <: Number, vType <: Function, fType <: Function,
            BType <: Function, GType <: Function, pType <: Union{Tuple,Nothing}, N} <: AbstractEquationPSDE{dType, tType}

    d::Int
    m::Int
    ni::Int
    ns::Int
    v::vType
    f::fType
    B::BType
    G::GType
    t₀::tType
    q₀::Array{dType, N}           #Initial condition: N=1 - single deterministic, N=2 - single random or multiple deterministic, N=3 - multiple deterministic
    p₀::Array{dType, N}
    parameters::pType
    periodicity::Vector{dType}

    function PSDE(m, ns, v::vType, f::fType, B::BType, G::GType,
                  t₀::tType, q₀::DenseArray{dType}, p₀::DenseArray{dType};
                  parameters=nothing, periodicity=zeros(dType,size(q₀,1))) where {
                        dType <: Number, tType <: Number,
                        vType <: Function, fType <: Function,
                        BType <: Function, GType <: Function}

        @assert size(q₀)  == size(p₀)
        @assert ndims(q₀) == ndims(p₀)

        N  = ndims(q₀)
        d  = size(q₀,1)
        ni = size(q₀,2)

        @assert N ∈ (1,2)
        @assert ni ≥ 1
        @assert ns ≥ 1
        @assert ni == 1 || ns == 1
        # either multiple deterministic initial conditions and one sample path
        # or one deterministic initial condition and multiple sample paths

        new{dType,tType,vType,fType,BType,GType,typeof(parameters),N}(d, m, ni, ns, v, f, B, G, t₀, q₀, p₀, parameters, periodicity)
    end
end

function PSDE(m::Int, ns::Int, v::Function, f::Function, B::Function, G::Function, q₀::DenseArray{DT}, p₀::DenseArray{DT}; kwargs...) where {DT}
    PSDE(m, ns, v, f, B, G, zero(DT), q₀, p₀; kwargs...)
end

Base.hash(sde::PSDE, h::UInt) = hash(sde.d, hash(sde.m, hash(sde.ni, hash(sde.ns,
                                hash(sde.v, hash(sde.f, hash(sde.B, hash(sde.G,
                                hash(sde.t₀, hash(sde.q₀, hash(sde.p₀,
                                hash(sde.parameters, hash(sde.periodicity, h)))))))))))))

Base.:(==)(sde1::PSDE, sde2::PSDE) = (
                                sde1.d == sde2.d
                             && sde1.m == sde2.m
                             && sde1.ni == sde2.ni
                             && sde1.ns == sde2.ns
                             && sde1.v == sde2.v
                             && sde1.f == sde2.f
                             && sde1.B == sde2.B
                             && sde1.G == sde2.G
                             && sde1.t₀ == sde2.t₀
                             && sde1.q₀ == sde2.q₀
                             && sde1.p₀ == sde2.p₀
                             && sde1.parameters  == sde2.parameters
                             && sde1.periodicity == sde2.periodicity)

function Base.similar(sde::PSDE, q₀::AbstractArray, p₀::AbstractArray, ns::Int)
    similar(sde, sde.t₀, q₀, p₀, ns)
end

function Base.similar(sde::PSDE, q₀::AbstractArray, p₀::AbstractArray)
    similar(sde, sde.t₀, q₀, p₀)
end

function Base.similar(sde::PSDE, t₀::TT, q₀::AbstractArray{DT,N}, p₀::AbstractArray{DT,N}, ns::Int=(N > 1 ? 1 : sde.ns)) where {DT <: Number, TT <: Number, N}
    @assert size(q₀) == size(p₀)
    @assert sde.d == size(q₀,1) == size(p₀,1)
    PSDE(sde.m, ns, sde.v, sde.f, sde.B, sde.G, t₀, q₀, p₀; parameters=sde.parameters, periodicity=sde.periodicity)
end

@inline Base.ndims(sde::PSDE) = sde.d

@inline periodicity(equation::PSDE) = equation.periodicity