text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Heat Equation

using NeuralPDE, Flux, ModelingToolkit, Optim, DiffEqFlux
using GalacticOptim
using Base
using SymPy

@parameters x t θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dt'~t

quasi_dirac = function(x,a)
    return 1/(a*pi^0.5)*exp(-(x/a)*(x/a))
end

quasi_dirac(0,0)

# 3D PDE
eq  = Dt(u(x,t,θ)) ~ Dxx(u(x,t,θ))

# Initial and boundary conditions
bcs = [u(x,0,θ) ~ quasi_dirac(x,1/100000)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,2.0),
           t ∈ IntervalDomain(0.0,2.0)]

# Discretization
dx = 0.25; dt = 0.25
# Neural network
chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = NeuralPDE.PhysicsInformedNN([dx,dt],
                                             chain,
                                             strategy = NeuralPDE.GridTraining())
pde_system = PDESystem(eq,bcs,domains,[x,t],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.01), progress = false; cb = cb, maxiters=1000)
phi = discretization.phi

ts,xs = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

u_real  = [[analytic_sol_func(t,x)[i] for t in ts for x in xs] for i in 1:3]
u_predict  = [[phi([t,x],res.minimizer)[i] for t in ts for x in xs] for i in 1:3]
diff_u = [abs.(u_real[i] .- u_predict[i] ) for i in 1:3]

for i in 1:3
    p1 = plot(xs, ts, u_real[i], st=:surface,title = "u$i, analytic");
    p2 = plot(xs, ts, u_predict[i], st=:surface,title = "predict");
    p3 = plot(xs, ts, diff_u[i],linetype=:contourf,title = "error");
    plot(p1,p2,p3)
    savefig("sol_u$i")
end


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

import Compat.Test
Test.@testset "Package ElasticArrays" begin
    include("util.jl")
    include("elasticarray.jl")
end

ext: jl


text:
module MultirateRungeKuttaMethod
using ..ODESolvers
using ..AdditiveRungeKuttaMethod
using ..LowStorageRungeKuttaMethod
using ..StrongStabilityPreservingRungeKuttaMethod
using ..MPIStateArrays: device, realview

using GPUifyLoops
include("MultirateRungeKuttaMethod_kernels.jl")

export MultirateRungeKutta

ODEs = ODESolvers
LSRK2N = LowStorageRungeKutta2N
SSPRK = StrongStabilityPreservingRungeKutta

"""
    MultirateRungeKutta(slow_solver, fast_solver; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f_fast(Q, t) + f_slow(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a multirate Runge-Kutta scheme using two different RK
solvers. This is based on

Currently only the low storage RK methods can be used as slow solvers

  - [`LowStorageRungeKuttaMethod`](@ref)

### References

    @article{SchlegelKnothArnoldWolke2012,
      title={Implementation of multirate time integration methods for air
             pollution modelling},
      author={Schlegel, M and Knoth, O and Arnold, M and Wolke, R},
      journal={Geoscientific Model Development},
      volume={5},
      number={6},
      pages={1395--1405},
      year={2012},
      publisher={Copernicus GmbH}
    }
"""
mutable struct MultirateRungeKutta{SS, FS, RT} <: ODEs.AbstractODESolver
  "slow solver"
  slow_solver::SS
  "fast solver"
  fast_solver::FS
  "time step"
  dt::RT
  "time"
  t::RT

  function MultirateRungeKutta(slow_solver::LSRK2N,
                               fast_solver,
                               Q=nothing;
                               dt=ODEs.getdt(slow_solver), t0=slow_solver.t
                              ) where {AT<:AbstractArray}
    SS = typeof(slow_solver)
    FS = typeof(fast_solver)
    RT = real(eltype(slow_solver.dQ))
    new{SS, FS, RT}(slow_solver, fast_solver, RT(dt), RT(t0))
  end
end
MRRK = MultirateRungeKutta

function MultirateRungeKutta(solvers::Tuple, Q=nothing;
                             dt=ODEs.getdt(solvers[1]), t0=solvers[1].t
                            ) where {AT<:AbstractArray}
  if length(solvers) < 2
    error("Must specify atleast two solvers")
  elseif length(solvers) == 2
    fast_solver = solvers[2]
  else
    fast_solver = MultirateRungeKutta(solvers[2:end], Q; dt = dt, t0=t0)
  end

  slow_solver = solvers[1]

  MultirateRungeKutta(slow_solver, fast_solver, Q; dt = dt, t0=t0)
end

function ODEs.dostep!(Q, mrrk::MultirateRungeKutta, param,
                      timeend::AbstractFloat, adjustfinalstep::Bool)
  time, dt = mrrk.t, mrrk.dt
  @assert dt > 0
  if adjustfinalstep && time + dt > timeend
    dt = timeend - time
    @assert dt > 0
  end

  ODEs.dostep!(Q, mrrk, param, time, dt)

  if dt == mrrk.dt
    mrrk.t += dt
  else
    mrrk.t = timeend
  end
  return mrrk.t
end

function ODEs.dostep!(Q, mrrk::MultirateRungeKutta{SS}, param,
                      time::AbstractFloat, dt::AbstractFloat,
                      in_slow_δ = nothing, in_slow_rv_dQ = nothing,
                      in_slow_scaling = nothing) where {SS <: LSRK2N}
  slow = mrrk.slow_solver
  fast = mrrk.fast_solver

  slow_rv_dQ = realview(slow.dQ)

  threads = 256
  blocks = div(length(realview(Q)) + threads - 1, threads)

  for slow_s = 1:length(slow.RKA)
    # Currnent slow state time
    slow_stage_time = time + slow.RKC[slow_s] * dt

    # Evaluate the slow mode
    slow.rhs!(slow.dQ, Q, param, slow_stage_time, increment = true)

    if in_slow_δ !== nothing
      slow_scaling = nothing
      if slow_s == length(slow.RKA)
        slow_scaling = in_slow_scaling
      end
      # update solution and scale RHS
      @launch(device(Q), threads=threads, blocks=blocks,
              update!(slow_rv_dQ, in_slow_rv_dQ, in_slow_δ, slow_scaling))
    end

    # Fractional time for slow stage
    if slow_s == length(slow.RKA)
      γ = 1 - slow.RKC[slow_s]
    else
      γ = slow.RKC[slow_s + 1] - slow.RKC[slow_s]
    end

    # RKB for the slow with fractional time factor remove (since full
    # integration of fast will result in scaling by γ)
    slow_δ = slow.RKB[slow_s] / (γ)

    # RKB for the slow with fractional time factor remove (since full
    # integration of fast will result in scaling by γ)
    nsubsteps = ODEs.getdt(fast) > 0 ? ceil(Int, γ * dt / ODEs.getdt(fast)) : 1
    fast_dt = γ * dt / nsubsteps

    for substep = 1:nsubsteps
      slow_rka = nothing
      if substep == nsubsteps
        slow_rka = slow.RKA[slow_s%length(slow.RKA) + 1]
      end
      fast_time = slow_stage_time + (substep - 1) * fast_dt
      ODEs.dostep!(Q, fast, param, fast_time, fast_dt, slow_δ, slow_rv_dQ,
                   slow_rka)
    end
  end
end

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractData{T,N}

Spatial data in a `N`-dimensional space with coordinates of type `T`.
"""
abstract type AbstractData{T,N} <: AbstractSpatialObject{T,N} end

"""
    variables(spatialdata)

Return the variable names in `spatialdata` and their types.
"""
variables(spatialdata::AbstractData) =
  Dict(var => eltype(array) for (var,array) in spatialdata.data)

"""
    spatialdata[ind,var]
    spatialdata[inds,vars]

Return the value of `var` for the `ind`-th point in `spatialdata`.
"""
Base.getindex(spatialdata::AbstractData, ind::Int, var::Symbol) =
  spatialdata.data[var][ind]

Base.getindex(spatialdata::AbstractData, inds::AbstractVector{Int}, var::Symbol) =
  [getindex(spatialdata, ind, var) for ind in inds]

Base.getindex(spatialdata::AbstractData, ind::Int, vars::AbstractVector{Symbol}) =
  [getindex(spatialdata, ind, var) for var in vars]

Base.getindex(spatialdata::AbstractData, inds::AbstractVector{Int}, vars::AbstractVector{Symbol}) =
  [getindex(spatialdata, ind, var) for ind in inds, var in vars]

"""
    spatialdata[var]

Return the values of `var` for all points in `spatialdata` with the shape
of the underlying domain.
"""
Base.getindex(spatialdata::AbstractData, var::Symbol) =
  [getindex(spatialdata, ind, var) for ind in 1:npoints(spatialdata)]

"""
    view(spatialdata, inds)
    view(spatialdata, vars)
    view(spatialdata, inds, vars)

Return a view of `spatialdata` with all points in `inds` and
all variables in `vars`.
"""
Base.view(sdata::AbstractData, inds::AbstractVector{Int}) =
  DataView(sdata, inds, collect(keys(variables(sdata))))

Base.view(sdata::AbstractData, vars::AbstractVector{Symbol}) =
  DataView(sdata, 1:npoints(sdata), vars)

Base.view(sdata::AbstractData, inds::AbstractVector{Int},
                               vars::AbstractVector{Symbol}) =
  DataView(sdata, inds, vars)

#----------------
# TABLES.JL API
#----------------
Tables.istable(::Type{<:AbstractData}) = true

Tables.columnaccess(::Type{<:AbstractData}) = true

function Tables.columns(spatialdata::AbstractData)
  vars = keys(variables(spatialdata))
  vals = [getindex(spatialdata, 1:npoints(spatialdata), var) for var in vars]
  NamedTuple{tuple(vars...)}(vals)
end

#-----------------
# MISSING VALUES
#-----------------
"""
    isvalid(spatialdata, ind, var)

Return `true` if the `ind`-th point in `spatialdata` has a valid value for `var`.
"""
function Base.isvalid(spatialdata::AbstractData, ind::Int, var::Symbol)
  val = getindex(spatialdata, ind, var)
  !(val ≡ missing || (val isa Number && isnan(val)))
end

"""
    valid(spatialdata, var)

Return all points in `spatialdata` with a valid value for `var`. The output
is a tuple with the matrix of coordinates as the first item and the vector
of values as the second item.
"""
function valid(spatialdata::AbstractData{T,N}, var::Symbol) where {N,T}
  V = variables(spatialdata)[var]
  npts = npoints(spatialdata)

  # pre-allocate memory for result
  X = Matrix{T}(undef, N, npts)
  z = Vector{V}(undef, npts)

  nvalid = 0
  for ind in 1:npoints(spatialdata)
    if isvalid(spatialdata, ind, var)
      nvalid += 1
      coordinates!(view(X,:,nvalid), spatialdata, ind)
      z[nvalid] = getindex(spatialdata, ind, var)
    end
  end

  X[:,1:nvalid], z[1:nvalid]
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, spatialdata::AbstractData{T,N}) where {N,T}
  npts = npoints(spatialdata)
  print(io, "$npts SpatialData{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", spatialdata::AbstractData{T,N}) where {N,T}
  println(io, spatialdata)
  println(io, "  variables")
  varlines = ["    └─$var ($V)" for (var,V) in variables(spatialdata)]
  print(io, join(sort(varlines), "\n"))
end

#------------------
# IMPLEMENTATIONS
#------------------
include("data/curve_data.jl")
include("data/geodataframe.jl")
include("data/point_set_data.jl")
include("data/regular_grid_data.jl")
include("data/structured_grid_data.jl")


text:
using StatsBase, DataFrames, Flux
include("2048_fn.jl")

using Flux
model = Flux.Chain(
  Dense(14*16+4, 128, relu),
  Dense(128, 64, relu),
  Dense(64, 32, relu),
  Dense(32, 1, relu));

loss(x,y) = sum((model(x) .- y).^2)
opt = ADAM(params(model))

@time simulate_one_game_and_train!(model, opt, 0.9)

init_grid = init_game()
@time (init_grid, fnl, seq, cart, tf, reward) = simulate_game!(init_grid, model, 0.5)

# score every move using untrained model
grid_try = copy(init_grid)
initial_model = [Tracker.data(model(grid_move_to_feature(sim_seq(grid_try, seq[1:i+1], cart[1:i+1], tf[1:i+1]), seq[i])))[1] for i=1:length(seq)-1]

# now fit the model once
x = generate_x_based_on_episode(init_grid, seq, cart, tf);
y = discounted_rewards(reward, 0.99)[1:end-1];
xy = [(x[:,j], y[j]) for j in 1:size(x)[2]]
Flux.@epochs 100 Flux.train!(loss, xy, opt)

updated_model = [Tracker.data(model(grid_move_to_feature(sim_seq(grid_try, seq[1:i+1], cart[1:i+1], tf[1:i+1]), seq[i])))[1] for i=1:length(seq)-1]
using Plots
plot(cumsum(reward))
plot!(updated_model)
plot!(y)


text:
@testset "hcat" begin
    @test isa(hcat(), TransposedVector{Any})
    @test isa(Base.typed_hcat(Int), TransposedVector{Int})

    @test isa([1 2 3], TransposedVector{Int})
    @test isa(Float64[1 2 3], TransposedVector{Float64})

    @test isa([1 2 3.0], TransposedVector{Float64})
    @test isa(Float32[1 2 3.0], TransposedVector{Float32})

    @test isa([[1 2 3] 4], TransposedVector{Int})
    @test isa([[1 2 3] [4 5]], TransposedVector{Int})

end


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

name = "PYTHIA"
version = v"8.244.0"

# Collection of sources required to complete build
sources = [
    "http://home.thep.lu.se/~torbjorn/pythia8/pythia8244.tgz" =>
    "e34880f999daf19cdd893a187123927ba77d1bf851e30f6ea9ec89591f4c92ca",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd pythia8244/
./configure --prefix=${prefix} --enable-shared --enable-64bit
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    Linux(:powerpc64le, libc=:glibc),
    MacOS(:x86_64)
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpythia8", :libpythia8)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
module PETScBinaryIO

export writepetsc, readpetsc

using SparseArrays

classids = Dict("Vec"=>1211214, "Mat"=>1211216, "IS"=>1211218)
ids_to_class = Dict(zip(values(classids), keys(classids)))

# PETSc IO for binary matrix IO. Format documentation here:
# http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatLoad.html
#  int    MAT_FILE_CLASSID
#  int    number of rows
#  int    number of columns
#  int    total number of nonzeros
#  int    *number nonzeros in each row
#  int    *column indices of all nonzeros (starting index is zero)
#  PetscScalar *values of all nonzeros

function write_(io :: IO, mat :: SparseMatrixCSC, int_type, scalar_type)
    m = SparseMatrixCSC(transpose(mat))
    rows, cols = size(mat)
    write(io, hton(int_type(classids["Mat"]))) # MAT_FILE_CLASSID
    write(io, hton(int_type(rows))) # number of rows
    write(io, hton(int_type(cols))) # number of columns
    write(io, hton(int_type(nnz(mat)))) # number of nonzeros

    # write row lengths
    for i = 1:rows
        write(io, hton(int_type(length(nzrange(m, i)))))
    end

    # write column indices
    cols = rowvals(m)
    for i = 1:rows
        for j in nzrange(m, i)
            write(io, hton(int_type(cols[j] - 1))) # PETSc uses 0-indexed arrays
        end
    end

    # write nonzero values
    vals = nonzeros(m)
    for i = 1:rows
        for j in nzrange(m, i)
            write(io, hton(scalar_type(vals[j])))
        end
    end
end

function write_(io :: IO, vec :: Vector, int_type, scalar_type)
    if eltype(vec) <: Integer
        write(io, hton(int_type(classids["IS"])))
    else
        write(io, hton(int_type(classids["Vec"])))
    end
    write(io, hton(int_type(length(vec)))) # number of rows
    if eltype(vec) <: Integer
        write(io, hton.(int_type.(vec .- 1)))
    else
        write(io, hton.(scalar_type.(vec)))
    end
end

"""
### writepetsc(filename, mat :: Vector{Union{SparseMatrixCSC, Vector}}; int_type = nothing, scalar_type = nothing)

Write a sparse matrix or vector to `filename` in a format PETSc can understand.
`PETSC_DIR` and `PETSC_ARCH` will be used to determine the size of types
written. `int_type` and `scalar_type` can be used to manually specify which
types are desired.
"""
function writepetsc(filename :: AbstractString, objs; int_type = nothing, scalar_type = nothing)
    int_type, scalar_type = determine_type(int_type, scalar_type)
    open(filename, "w") do io
        for o in objs
            write_(io, o, int_type, scalar_type)
        end
    end
end

function writepetsc(filename :: AbstractString, mat :: SparseMatrixCSC; int_type = nothing, scalar_type = nothing)
    writepetsc(filename, [mat]; int_type=int_type, scalar_type=scalar_type)
end

function writepetsc(filename :: AbstractString, vec :: Vector{T}; int_type = nothing, scalar_type = nothing) where T <: Number
    writepetsc(filename, [vec]; int_type=int_type, scalar_type=scalar_type)
end

function read_prefix_vec(io, int_type, scalar_type)
    len = ntoh(read(io, int_type))
    read_vec(io, scalar_type, len)
end

function read_vec(io, ty, sz)
    ary = Array{ty}(undef, sz)
    read!(io, ary)
    ntoh.(ary)
end

function read_mat(io, int_type, scalar_type)
    rows = ntoh(read(io, int_type))
    cols = ntoh(read(io, int_type))
    nnz = ntoh(read(io, int_type))

    row_ptr = Array{int_type}(undef, rows+1)
    row_ptr[1] = 1

    # read row lengths
    row_ptr[2:end] = read_vec(io, int_type, rows)
    cumsum!(row_ptr, row_ptr)

    # write column indices
    colvals = read_vec(io, int_type, nnz) .+ int_type(1)

    # write nonzero values
    vals = read_vec(io, scalar_type, nnz)

    mat = SparseMatrixCSC(cols, rows, row_ptr, colvals, vals)
    SparseMatrixCSC(transpose(mat))
end

function read_single(io, int_type, scalar_type)
    class_id = ntoh(read(io, int_type))
    if !in(class_id, keys(ids_to_class))
        throw("Invalid PETSc binary file $class_id")
    end
    if ids_to_class[class_id] == "Vec"
        read_prefix_vec(io, int_type, scalar_type)
    elseif ids_to_class[class_id] == "Mat"
        read_mat(io, int_type, scalar_type)
    elseif ids_to_class[class_id] == "IS"
        read_prefix_vec(io, int_type, int_type) .+ 1
    else
        error("Invalid class id $class_id")
    end
end

"""
Tries to find and parse a petscvariables file.
"""
function parse_petsc_config()
    default_sizes = (Int32, Float64)
    isizemap = Dict(32 => Int32, 64 => Int64)
    ssizemap = Dict(32 => Float32, 64 => Float64)
    path = nothing
    if haskey(ENV, "PETSC_DIR")
        dir = ENV["PETSC_DIR"]
        path = joinpath(dir, "lib", "petsc", "conf", "petscvariables")
        if haskey(ENV, "PETSC_ARCH")
            arch = ENV["PETSC_ARCH"]
            path = joinpath(dir, arch, "lib", "petsc", "conf", "petscvariables")
        end
    end

    if path != nothing && isfile(path)
        isize = -1
        ssize = -1
        for line in eachline(path)
            if occursin("PETSC_SCALAR_SIZE", line)
                ssize = parse(Int, match(r"[0-9]+", line).match)
            end
            if occursin("PETSC_INDEX_SIZE", line)
                isize = parse(Int, match(r"[0-9]+", line).match)
            end
        end
        if !haskey(isizemap, isize) || !haskey(ssizemap, ssize)
            @warn "Could not determine PETSC_INDEX_SIZE and PETSC_SCALAR_SIZE from $path, defaulting to $default_sizes" maxlog=1
            return default_sizes
        end
        return (isizemap[isize], ssizemap[ssize])
    elseif path != nothing
        @warn "$path does not exist, defaulting to $default_sizes" maxlog=1
    end
    default_sizes
end

function determine_type(int_type, scalar_type)
    config_int_type, config_scalar_type = parse_petsc_config()
    if int_type == nothing
        int_type = config_int_type
    end
    if scalar_type == nothing
        scalar_type = config_scalar_type
    end
    int_type, scalar_type
end

"""
### readpetsc(filename; int_type = nothing, scalar_type = nothing)
           :: Vector{Union{SparseMatrixCSC, Vector}}

Read a sparse matrix in PETSc's binary format from `filename`. PETSC_DIR and
PETSC_ARCH will be used to determine the correct integer and float type to use.
`int_type` and `scalar_type` can be used to override these values.
"""
function readpetsc(filename; int_type = nothing, scalar_type = nothing) :: Vector{Union{SparseMatrixCSC, Vector}}
    int_type, scalar_type = determine_type(int_type, scalar_type)
    open(filename) do io
        items = []
        while !eof(io)
            push!(items, read_single(io, int_type, scalar_type))
        end
        items
    end
end

end


text:
module SphericalFunctions

using FastTransforms, LinearAlgebra, ProgressMeter
import Base.Threads: @threads, nthreads

export complex_powers, complex_powers!
export theta_phi, phi_theta
export WignerMatrixCalculator, H!, d!, D!, Y!
export ALFRecursionCoefficients, ALFrecurse!, ALFcompute!
export fejer1, fejer2, clenshaw_curtis
export map2salm, map2salm!, plan_map2salm

const MachineFloat = Union{Float16, Float32, Float64}

include("utils.jl")

include("complex_powers.jl")

include("indexing.jl")

include("wigner_matrices/evaluate.jl")

include("associated_legendre/calculator.jl")
using .AssociatedLegendreFunction: ALFRecursionCoefficients, ALFrecurse!, ALFcompute!

include("weights.jl")

include("map2salm.jl")

end # module


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    PredicatePartitioner(pred)

A method for partitioning spatial objects with a given predicate
function `pred(i, j)`.
"""
struct PredicatePartitioner <: AbstractPredicatePartitioner
  pred::Function
end

(p::PredicatePartitioner)(i, j) = p.pred(i, j)


text:
module Canvas
using Printf
using Markdown
using Cairo
using ComplexVisual
@ComplexVisual.import_huge

import Main.DocGenerator: DocSource, DocCreationEnvironment, DocContext,
        Document, substitute_marker_in_markdown, create_doc_icon, append_md

"""
# [![./Canvas_docicon.png]({image_from_canvas: get_doc_icon()}) Canvas](./Canvas.md)

Canvases (i.e. subtypes of `CV_Canvas`) are thin wrappers for Cairo canvases.

## `doc: CV_Canvas`

## `doc: CV_2DCanvas`

## `doc: CV_Std2DCanvas`

## `doc: CV_Math2DCanvas`

## `doc: cv_math2pixel:Tuple{CV_Math2DCanvas, Float64, Float64}`

## `doc: cv_pixel2math:Tuple{CV_Math2DCanvas, Integer, Integer}`

## `doc: cv_save_image`
"""
canvas_intro() = nothing

function get_doc_icon()
    layout = CV_2DLayout()

    can = cv_filled_canvas(30, 40, cv_white)
    can_l = cv_add_canvas!(layout, can, cv_anchor(can, :center), (0, 0))

    r_rect = cv_filled_canvas(18, 3, cv_color(1, 0, 0, 0.7))
    r_rect_l = cv_add_canvas!(layout, r_rect,
        cv_anchor(r_rect, :west),
        cv_translate(cv_anchor(can_l, :west), 3, -4))

    b_rect = cv_filled_canvas(12, 18, cv_color(0, 0, 1, 0.7))
    b_rect_l = cv_add_canvas!(layout, b_rect,
        cv_anchor(b_rect, :center),
        cv_translate(cv_anchor(can_l, :center), 7, 0))

    bg_fill = (cv_color(0, 0, 0, 0) → cv_op_over) ↦ CV_FillPainter()
    stand_canvas = CV_Math2DCanvas(-1.0 + 1.0im, 1.0 - 0.7im, 10)
    line_segs = [
        [-0.2 + 1.0im, -0.7 - 0.7im],
        [+0.2 + 1.0im, +0.7 - 0.7im],
        [0.0 + 1.0im,  0.0 - 0.7im]]

    lp_style = cv_color(0.51, 0.16, 0) → cv_linewidth(2)
    lp = lp_style ↦ CV_LinePainter(line_segs)

    cv_create_context(stand_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, lp)
    end

    stand_canvas_l = cv_add_canvas!(layout, stand_canvas,
        cv_anchor(stand_canvas, :north), cv_anchor(can_l, :south))

    stand_border = cv_border(layout, can_l, 0, 0, 2, 0; style=lp_style)
    shad_border = cv_border(layout, can_l, 2, 2, 0, 0;
        style=cv_color(0.4, 0.4, 0.4))

    cv_ensure_size!(layout, 70, 70)

    can_layout = cv_canvas_for_layout(layout)
    con_layout = cv_create_context(can_layout;
        fill_with=cv_color(0.7, 0.7, 0.7))
    can_l(con_layout)
    stand_canvas_l(con_layout)
    stand_border(con_layout)
    shad_border(con_layout)
    r_rect_l(con_layout)
    b_rect_l(con_layout)
    cv_destroy(con_layout)

    icon = create_doc_icon(can_layout)
    return icon
end

function create_document(doc_env::DocCreationEnvironment)
    doc_source = DocSource("Canvas", @__MODULE__)
    context = DocContext(doc_env, doc_source)

    md = Markdown.MD()
    for part in (canvas_intro, )
        part_md = Base.Docs.doc(part)
        substitute_marker_in_markdown(context, part_md)

        append_md(md, part_md)
    end

    doc = Document(doc_source, md)
    return doc
end

end
# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
module OVERTVerify

# include stuff
include("models/problems.jl")
include("overt_to_mip.jl")

## The following code in the nv folder forked from https://github.com/sisl/NeuralVerification.jl
include("nv/utils/activation.jl")
include("nv/utils/network.jl")
include("nv/utils/problem.jl")
include("nv/utils/util.jl")
include("nv/optimization/utils/constraints.jl")
include("nv/optimization/utils/objectives.jl")
include("nv/optimization/utils/variables.jl")
include("nv/reachability/maxSens.jl")
##

include("logic.jl")
include("reachability_functions.jl")
include("models/single_pendulum/single_pend.jl")
include("models/tora/tora.jl")
include("models/car/simple_car.jl")
include("models/acc/acc.jl")
include("overt_parser_minimal.jl")
include("dreal_wrapper/smt2_interface.jl")
include("dreal_wrapper/proof_functions.jl")
include("dreal_wrapper/dreal_utils.jl")
include("dreal_wrapper/compare_to_dreal.jl")

export OvertQuery,
    SinglePendulum,
    ACC,
    SimpleCar,
    Tora,
    symbolic_satisfiability,
    symbolic_reachability_with_concretization,
    clean_up_sets,
    clean_up_meas_sets,
    check_avoid_set_intersection,
    compare_to_dreal,
    Id,
    ReLU,
    Constraint
        # stuff
end


text:
"""
Julia program to check if a number is a Strong number or not.
Strong Number is a number whose sum of the factorial of digits
is equal to the number itself.
"""

# Helper function to return the factorial of a number
function factorial(num)
    if ( num == 1 || num == 0 )
        return 1
    end
    return num * factorial( num - 1)
end

""" Extract each digit of the given number and calculate the sum of the
    factorial of each digit present in the given number"""
function strong_number(num)
    temp = num
    sum = 0
    while(temp > 0)
        rem = temp % 10
        sum = sum + factorial(rem)
        temp = temp ÷ 10
    end
    return sum
end

print("Enter the number: ")
num = readline()
num = parse(Int, num)
sum = strong_number(num)
if sum == num
    println("The given number $num is a Strong Number.")
else
    println("The given number $num is not a Strong Number.")
end

"""
Time Complexity: O(n*log(n)), where 'n' is the given number
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE 1
Enter the number: 3124
The given number 3124 is not a Strong Number.

SAMPLE 2
Enter the number: 145
The given number 145 is a Strong Number.
"""


text:
using Roots #use: pkg.add("Roots") or in package mode: add Roots
rhoW = 1000 #Density of water in kg/m^3
rhoPVC = 1380 #Density of PVC in kg/m^3
Cd = 0.7 #Drag coefficient
Vs = (4/3)*pi #Coefficient for volume of sphere
g = 9.81 #Acceleration by gravity on Earth
V_dot = 3.2e-3 #Volumetric flow rate
A_in = (pi/4)*0.15^2 #Area of inlet
A_out = (pi/4)*0.35^2 #Area of outlet
f(r) = (sqrt((2*(rhoPVC-rhoW)*g*(Vs*r^3))/
(Cd*rhoW*(pi*r^2))))*(A_in-(pi*r^2))-V_dot
    #assumes ball diameter is closer to inlet
f2(r)= (sqrt((2*(rhoPVC-rhoW)*g*(Vs*r^3))/
(Cd*rhoW*(pi*r^2))))*(A_out-(pi*r^2))-V_dot
    #assumes ball diameter is closer to outlet
a = 2*fzero(f,0.07) #reference Roots in Julia documentation for usage of fzero
display(a) #lower bound result
b = 2*fzero(f2,0.17)
display(b) #upper bound result


text:
using NodeCall

cd(@__DIR__)
# Install `node-fetch` with `npm`.
run(npm("install"))

const jsfetch = require("node-fetch")
const console = node"console"

# `fetch` is Julia's `await`
response = fetch(jsfetch("https://ip.sunoru.com"))

# Or you can use the await macro
ip = @await response.text()

@show ip

# Define an asynchronous function in JS scripts.
# Note that the imported `fetch` should be manually passed in.
f1 = node"""async (fetch) => {
    const response = await fetch('https://httpbin.org/post', {method: 'POST', body: 'a=1'})
    const data = await response.json()
    return data
}"""

p = f1(jsfetch).then() do data
    # We can use Julia's do-block to pass the function to `then()`
    @show data
    console.log(data)
    for k in keys(data)
        v = data[k]
        if v isa JsObject
            v = value(Dict, v)
        end
        @show k, v
    end
    throw("Some error in promise")
end.catch() do err
    # The `catch()` method is similar.
    @show err
    console.log("Catched: \n", err)
end

# Should explicitly wait for the promise to complete.
wait(p)


text:
"""
	TimLaneChanger
A simple lane change behavior that changes lanes whenever the lead car is going slower than our desired speed.
Lane changes are made when there is an available lane, fore/rear gaps exceed our thresholds, we are faster
than a rear vehicle in the target lane, and any lead vehicle in the target lane is faster than we can currently go.

Has not been published anywhere, so first use in a paper would have to describe this.
See MOBIL if you want a lane changer you can cite.

# Constructors
	TimLaneChanger(timestep::Float64;v_des::Float64=29.0,rec::SceneRecord=SceneRecord(2,timestep),threshold_fore::Float64 = 50.0,threshold_lane_change_gap_fore::Float64 = 10.0, threshold_lane_change_gap_rear::Float64 = 10.0,dir::Int=DIR_MIDDLE)

# Fields
- `dir::Int = DIR_MIDDLE` the desired lane to go to eg: left,middle (i.e. stay in same lane) or right
- `rec::SceneRecord` TODO
- `v_des::Float64 = 29.0` desired velocity
- `threshold_fore::Float64 = 50.0` Distance from lead vehicle
- `threshold_lane_change_gap_fore::Float64 = 10.0` Space in front
- `threshold_lane_change_gap_rear::Float64 = 10.0` Space rear
"""
mutable struct TimLaneChanger <: LaneChangeModel
    dir::Int
    rec::SceneRecord

    v_des::Float64
    threshold_fore::Float64
    threshold_lane_change_gap_fore::Float64
    threshold_lane_change_gap_rear::Float64

    function TimLaneChanger(
        timestep::Float64;
        v_des::Float64=29.0,
        rec::SceneRecord=SceneRecord(2,timestep),
        threshold_fore::Float64 = 50.0,
        threshold_lane_change_gap_fore::Float64 = 10.0,
        threshold_lane_change_gap_rear::Float64 = 10.0,
        dir::Int=DIR_MIDDLE,
        )

        retval = new()
        retval.dir = dir
        retval.rec = rec
        retval.v_des = v_des
        retval.threshold_fore = threshold_fore
        retval.threshold_lane_change_gap_fore = threshold_lane_change_gap_fore
        retval.threshold_lane_change_gap_rear = threshold_lane_change_gap_rear
        retval
    end
end
get_name(::TimLaneChanger) = "TimLaneChanger"
function set_desired_speed!(model::TimLaneChanger, v_des::Float64)
    model.v_des = v_des
    model
end
function observe!(model::TimLaneChanger, scene::Scene, roadway::Roadway, egoid::Int)

    rec = model.rec
    update!(rec, scene)
    vehicle_index = findfirst(egoid, scene)

    veh_ego = scene[vehicle_index]
    v = veh_ego.state.v

    left_lane_exists = convert(Float64, get(N_LANE_LEFT, rec, roadway, vehicle_index)) > 0
    right_lane_exists = convert(Float64, get(N_LANE_RIGHT, rec, roadway, vehicle_index)) > 0
    fore_M = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront(), max_distance_fore=model.threshold_fore)
    fore_L = get_neighbor_fore_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointRear(), VehicleTargetPointFront())
    fore_R = get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointRear(), VehicleTargetPointFront())
    rear_L = get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())
    rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

    model.dir = DIR_MIDDLE
    if fore_M.Δs < model.threshold_fore # there is a lead vehicle
        veh_M = scene[fore_M.ind]
        speed_M = veh_M.state.v
        if speed_M ≤ min(model.v_des, v) # they are driving slower than we want

            speed_ahead = speed_M

            # consider changing to a different lane
            if right_lane_exists &&
               fore_R.Δs > model.threshold_lane_change_gap_rear && # there is space rear
               rear_R.Δs > model.threshold_lane_change_gap_fore && # there is space fore
               (rear_R.ind == nothing || scene[rear_R.ind].state.v ≤ v) && # we are faster than any follower
               (fore_R.ind == nothing || scene[fore_R.ind].state.v > speed_ahead) # lead is faster than current speed

                speed_ahead = fore_R.ind != nothing ? scene[fore_R.ind].state.v : Inf
                model.dir = DIR_RIGHT
            end
            if left_lane_exists &&
               fore_L.Δs > model.threshold_lane_change_gap_rear && # there is space rear
               rear_L.Δs > model.threshold_lane_change_gap_fore && # there is space fore
               (rear_L.ind == nothing || scene[rear_L.ind].state.v ≤ v) && # we are faster than any follower
               (fore_L.ind == nothing || scene[fore_L.ind].state.v > speed_ahead) # lead is faster than current speed

                speed_ahead = fore_L.ind != nothing ? scene[fore_L.ind].state.v : Inf
                model.dir = DIR_LEFT
            end
        end
    end

    model
end
Base.rand(model::TimLaneChanger) = LaneChangeChoice(model.dir)


text:
#!/usr/bin/env julia

# Author: Hiroharu Sugawara <hsugawa@gmail.com>
# Copyright: (C) 2020 Hiroharu Sugawara

using PandocFiltersLiveJuliaCode

if abspath(PROGRAM_FILE) == @__FILE__
   PandocFiltersLiveJuliaCode.filter([
      PandocFiltersLiveJuliaCode.investigate_image,
      PandocFiltersLiveJuliaCode.investigate_codeblock ])
end


text:
__precompile__(true)
module CSV

using DataStreams, DataFrames, WeakRefStrings

export Data, DataFrame

if !isdefined(Core, :String)
    typealias String UTF8String
end

if Base.VERSION < v"0.5.0-dev+4631"
    unsafe_wrap{A<:Array}(::Type{A}, ptr, len) = pointer_to_array(ptr, len)
    unsafe_string(ptr, len) = utf8(ptr, len)
    unsafe_wrap(::Type{String}, ptr, len) = unsafe_string(ptr, len)
    escape_string(io, str1, str2) = print_escaped(io, str1, str2)
end

immutable CSVError <: Exception
    msg::String
end

const RETURN  = UInt8('\r')
const NEWLINE = UInt8('\n')
const COMMA   = UInt8(',')
const QUOTE   = UInt8('"')
const ESCAPE  = UInt8('\\')
const PERIOD  = UInt8('.')
const SPACE   = UInt8(' ')
const TAB     = UInt8('\t')
const MINUS   = UInt8('-')
const PLUS    = UInt8('+')
const NEG_ONE = UInt8('0')-UInt8(1)
const ZERO    = UInt8('0')
const TEN     = UInt8('9')+UInt8(1)
Base.isascii(c::UInt8) = c < 0x80

@inline function unsafe_read(from::Base.AbstractIOBuffer, ::Type{UInt8}=UInt8)
    @inbounds byte = from.data[from.ptr]
    from.ptr = from.ptr + 1
    return byte
end
unsafe_read(from::IO, T) = Base.read(from, T)

@inline function unsafe_peek(from::Base.AbstractIOBuffer)
    @inbounds byte = from.data[from.ptr]
    return byte
end
unsafe_peek(from::IO) = (mark(from); v = Base.read(from, UInt8); reset(from); return v)

"""
Represents the various configuration settings for csv file parsing.

Keyword Arguments:

 * `delim`::Union{Char,UInt8} = how fields in the file are delimited
 * `quotechar`::Union{Char,UInt8} = the character that indicates a quoted field that may contain the `delim` or newlines
 * `escapechar`::Union{Char,UInt8} = the character that escapes a `quotechar` in a quoted field
 * `null`::String = indicates how NULL values are represented in the dataset
 * `dateformat`::Union{AbstractString,Dates.DateFormat} = how dates/datetimes are represented in the dataset
"""
type Options
    delim::UInt8
    quotechar::UInt8
    escapechar::UInt8
    separator::UInt8
    decimal::UInt8
    null::String
    nullcheck::Bool
    dateformat::Dates.DateFormat
    datecheck::Bool
    # non-public for now
    datarow::Int
    rows::Int
    header::Union{Integer,UnitRange{Int},Vector}
    types::Union{Dict{Int,DataType},Dict{String,DataType},Vector{DataType}}
end

Options(;delim=COMMA,quotechar=QUOTE,escapechar=ESCAPE,null=String(""),dateformat=Dates.ISODateFormat, datarow=-1, rows=0, header=1, types=DataType[]) =
    Options(delim%UInt8,quotechar%UInt8,escapechar%UInt8,COMMA,PERIOD,
            null,null != "",isa(dateformat,Dates.DateFormat) ? dateformat : Dates.DateFormat(dateformat),dateformat == Dates.ISODateFormat, datarow, rows, header, types)
function Base.show(io::IO,op::Options)
    println(io,"    CSV.Options:")
    println(io,"        delim: '",Char(op.delim),"'")
    println(io,"        quotechar: '",Char(op.quotechar),"'")
    print(io,"        escapechar: '"); escape_string(io,string(Char(op.escapechar)),"\\"); println(io,"'")
    print(io,"        null: \""); escape_string(io,op.null,"\\"); println(io,"\"")
    print(io,"        dateformat: ",op.dateformat)
end

"""
constructs a `CSV.Source` file ready to start parsing data from

implements the `Data.Source` interface for providing convenient `Data.stream!` methods for various `Data.Sink` types
"""
type Source{I<:IO} <: Data.Source
    schema::Data.Schema
    options::Options
    data::I
    datapos::Int # the position in the IOBuffer where the rows of data begins
    fullpath::String
end

function Base.show(io::IO,f::Source)
    println(io,"CSV.Source: ",f.fullpath)
    println(io,f.options)
    show(io, f.schema)
end

"""
constructs a `CSV.Sink` file ready to start writing data to

implements the `Data.Sink` interface for providing convenient `Data.stream!` methods for various `Data.Source` types
"""
type Sink <: Data.Sink
    schema::Data.Schema
    options::Options
    path::String
    data::IO
    datapos::Int # the byte position in `io` where the data rows start
    header::Bool
    # quotefields::Bool # whether to always quote string fields or not
end

include("parsefields.jl")
include("io.jl")
include("Source.jl")
include("Sink.jl")

end # module