text:
function adaptive_sampling(A,b,obj,noisyp,Atest,btest,testx,instance,S)
    Experiment = 1

    I = 100
	TS = 100
    re_solve_counter = 0
    cone = fill(zeros(n,n), 1)
    Pred_Error_Evolution = zeros(Float64, I, 10)
    ind = 0

    A_train = zeros(Float64, 2*n+1, n, I)
    b_train = zeros(Float64, I, 2*n+1)
    x_train = zeros(Float64, 30, n, I)

    # Generate random folders to store BARON temp files
    id=randstring(5)
    @eval @everywhere id=$id
    mkdir(string(root,"Results//ex",id))
    for t=1:S
        mkdir(string(root,"Results//ex",id,"//opts",t))
    end

    while(Experiment <= I)
        T = Experiment
        A_train[:,:,Experiment] = deepcopy(A);
        # Scale A and b to use a smaller big-M parameter while solving (P1)
        A_train[n+1, :, iteration] = A_train[n+1, :, Experiment]/10000;
        b_train[iteration, :] = deepcopy(b);
        b_train[iteration, n+1] = b_train[iteration, n+1]/10000;
        x = SolveModel(A, obj, b, n)
        J = 30;
        noisyx = zeros(Float64, J, n)
        for j = 1:J
            for k = 1:n
                noisyx[j, k] = x[k] + randn(rng, Float64)/100
            end
        end
        x_train[:,:,iteration] = deepcopy(noisyx);

        hat_x = zeros(Float64, I, n)
        hat_z = zeros(Float64, I, c)

        hat_x[T, :], r, z = SolveP1(A_train[:,:,T], b_train[T,:], noisyx, 1, J)
        hat_z[T, :] = z
        istight = []
        for i = 1:c
            if isapprox(z[1, i], 1, atol = 1e-3) == 1
                push!(istight, i)
            end
        end

        if Experiment == 1
            for j = 1:length(istight)
                cone[iteration][j,:] = A[istight[j], :]/norm(A[istight[j], :])
            end
        else
            B = zeros(Float64, length(istight), n)
            for j = 1:length(istight)
                B[j,:] = A[istight[j], :]/norm(A[istight[j], :])
            end
            push!(cone, B)
        end
        status = SolveFP(cone, T)

        if status == :Infeasible
            re_solve_counter += 1
            z = []
            x, r, z, status = SolveP1m(A_train[:,:,1:T], b_train[1:T,:], x_train[:,:,1:T], T, J, instance, hat_x, hat_z)

            if status != :Optimal
                break
            end
            hat_x[1:T, :] = x
            hat_z[1:T, :] = z

            cone = fill(zeros(n,n), 1)

            for k = 1:Experiment
                istight = []
                for i = 1:c
                    if isapprox(z[k, i], 1, atol = 1e-3)
                        push!(istight, i)
                    end
                end
                if k == 1
                    for j = 1:length(istight)
                        cone[k][j,:] = A_train[istight[j], :, k]/norm(A_train[istight[j], :, k])
                    end
                else
                    B = zeros(Float64, length(istight), n)
                    for j = 1:length(istight)
                        B[j,:] = A_train[istight[j], :, k]/norm(A_train[istight[j], :, k])
                    end
                    push!(cone, B)
                end
            end
        end

        hatc1, loss1 = SolveP2(noisyp, cone, iteration)
        xhat  = zeros(Float32, TS, n)
        for test = 1:TS
            xhat[test, 1:n] = SolveModel(Atest[:, :, test], hatc1, btest[test, :], n)
        end
        errorrate = accuracy(testx, xhat)
        push!(Pred_Error_Evolution,errorrate)

        # Initialize arrays to store samples for heuristic-based adaptive sampling
        global bprime = zeros(Float64, S, c)
        global Aprime = zeros(Float64, c, n, S)
        xprime = zeros(Float64, S, dim)
        global zprime = zeros(Int16, S, c)
        phi = zeros(Float64, S)  # To store the objective values of (ASP)
        Boun = zeros(Float64, S) # To store the best bound for (ASP)
        G = zeros(Float64, S)

        for T = 1:S
            bprime[T, :] = b
            Aprime[:, :, T] = A
            for i = 1:dim
                Aprime[dim+1, i, T] = rand(rng, 5000:15000)/100
            end
            hatc2 = SampleC(cone) # Sample a random vector from the admissible set
            xprime[T, :] = SolveModel(Aprime[:,:,T], hatc2, bprime[T,:], n)
            temp, zprime[T, :] = tightind(Aprime[:,:,T], bprime[T,:], xprime[T,:])
        end

        # Solve the adaptive sampling problem for the samples in Aprime, bprime
        Results = zeros(Float64, S, 3)
        Results = pmap(T->maxmin(T, cone, Aprime[:,:,T], zprime[T,:]), 1:S, retry_delays=ones(10))
        for i = 1:S
            phi[i] = Results[i][1]
            Boun[i] = Results[i][2]
        end
        ind = findmax(phi)

        # Update A to the sample with the highest value of \eta
        A = Aprime[:,:,ind[2]]
        Experiment += 1
    end

    cd(string(root,"Results"))
    df1 = DataFrame(Predic_Error = Pred_Error_Evolution)
    s1 = @sprintf "Adaptive_Parallel_dual_%03.d_%03.d_%04.d.csv" instance dim S
    CSV.write(s1, df1)
end

max_stars_repo_name: ddolab/InvLP

max_issues_repo_name: ddolab/InvLP

max_forks_repo_name: ddolab/InvLP


text:
module TestPkg

1==1

end

ext: jl


text:
@testset "1d derivative" begin
  xs = collect(0:0.1:10)
  fs = [sin(x) for x in xs]

  m = mwls_kd(xs, fs, 1, (d, e) -> exp(-d^2))
  cll = mwls_cll(xs, fs, 1, (d, e) -> exp(-d^2))

  teststatus = true
  for p in 0:0.1:10
    isapprox(mwls_diff(m, p, 1), cos(p), atol=0.1)
  end
  @test teststatus == true

  teststatus = true
  for p in 0:0.1:10
    teststatus &= isapprox(mwls_diff(cll, p, 1), cos(p), rtol=0.1)
  end
  @test teststatus == true

  # not enough data for approximation, zero returned by default
  @test isapprox(mwls_diff(m, -0.9, 1), 0, atol=0.1)
  @test isapprox(mwls_diff(cll, -0.9, 1), 0, atol=0.1)
  @test isapprox(mwls_diff(m, 10.9, 1), 0, atol=0.1)
  @test isapprox(mwls_diff(cll, 10.9, 1), 0, atol=0.1)
end


text:
module Laplace

using StaticArrays
using LinearAlgebra
using ..AeroMDAO: Point2D, Point3D

abstract type AbstractLaplace end

# Performs velocity and potential computations for an object on a grid
grid_data(object :: AbstractLaplace, xs) = velocity(object, xs), potential(object, xs)
velocity(object :: AbstractLaplace, xs) = map(x -> velocity(object, x...), xs)
potential(object :: AbstractLaplace, xs) = map(x -> potential(object, x...), xs)

# Performs velocity and potential calculations on a grid
function grid_data(objects :: Vector{<: AbstractLaplace}, xs)
    vels = foldl((v1, v2) -> [ u .+ v for (u, v) ∈ zip(v1, v2) ], velocity(object, xs) for object ∈ objects)
    pots = foldl((v1, v2) -> v1 + v2, potential(object, xs) for object ∈ objects)

    vels, pots
end

## 2D singularities
#============================================#

struct Singularity2D{T <: Real} <: AbstractLaplace
    strength :: T
    r 		 :: Point2D{T}
end

# Getters
strength(s :: Singularity2D) = s.strength
x(s :: Singularity2D) 		 = s.r.x
y(s :: Singularity2D) 		 = s.r.y

source_velocity(src :: Singularity2D, x, y) = SVector(strength(src) / (2π) * (x - y(src)) / ((x - y(src))^2 + (y - x(src))^2), str / (2π) * (y - x(src)) / ((x - y(src))^2 + (y - x(src))^2))
source_potential(src :: Singularity2D, x, y) = strength(src) / (4π) * log((x - y(src))^2 + (y - x(src))^2)
source_stream(src :: Singularity2D, x, y) = strength(src) / (2π) * atan(y - x(src), x - y(src))

doublet_velocity(dub :: Singularity2D, x, y) = SVector(strength(dub) / (2π) * ((x - y(dub))^2 - (y - x(dub))^2) / ((x - y(dub))^2 + (y - x(dub))^2)^2, - strength(dub) / (2π) * 2 * (x - y(dub)) * (y - x(dub)) / ((x - y(dub))^2 + (y - x(dub))^2)^2)
doublet_potential(dub :: Singularity2D, x, y) = -strength(dub) / (2π) * (y - x(dub)) / ((x - y(dub))^2 + (y - x(dub))^2)
doublet_stream(dub :: Singularity2D, x, y) = -strength(dub) / (2π) * (y - x(dub)) / ((x - y(dub))^2 + (y - x(dub))^2)

vortex_velocity(vor :: Singularity2D, x, y) = SVector(-strength(vor) / (2π) * (y - x(vor)) / ((x - y(vor))^2 + (y - x(vor))^2), str / (2π) * (x - y(vor)) / ((x - y(vor))^2 + (y - x(vor))^2))
vortex_potential(vor :: Singularity2D, x, y) = strength(vor) / (2π) * atan(y - x(vor), x - y(vor))
vortex_stream(vor :: Singularity2D, x, y) = -strength(vor) / (4π) * log((x - y(vor))^2 + (y - x(vor))^2)

struct Uniform2D{T <: Real} <: AbstractLaplace
    magnitude :: T
    angle 	  :: T
    Uniform2D{T}(mag, ang) where T <: Real = new(mag, deg2rad(ang))
end

magnitude(uni :: Uniform2D) = uni.magnitude
angle(uni :: Uniform2D)     = uni.angle

Uniform2D(mag :: T, ang :: T) where T <: Real = Uniform2D{T}(mag, ang)
Uniform2D(mag, ang) = Uniform2D(promote(mag, ang)...)

velocity(uni :: Uniform2D) = let (sa, ca) = sincos(uni.angle); uni.magnitude * SVector(ca, sa) end
potential(uni :: Uniform2D, x, y) = uni.magnitude * (x * cos(uni.angle) + y * sin(uni.angle))
stream(uni :: Uniform2D, x, y)    = uni.magnitude * (y * cos(uni.angle) - x * sin(uni.angle))

## 3D singularities
#============================================#

# struct Singularity3D{T <: Real} <: AbstractLaplace
# 	str :: T
# 	r 	:: Point3D{T}
# end

# source_velocity(src :: Source2D, x, y, z)
# source_potential(src :: Source2D, x, y, z)
# source_stream(src :: Source2D, x, y, z)

## Freestream
#============================================#

struct Freestream{T <: Real} <: AbstractLaplace
    V 	  :: T
    alpha :: T
    beta  :: T
    omega :: SVector{3,T}
    Freestream{T}(V, α_deg, β_deg, Ω) where T <: Real = new(V, deg2rad(α_deg), deg2rad(β_deg), Ω)
end

"""
    Freestream(V, α, β, Ω)

A Freestream flow in spherical polar coordinates with magnitude ``V``, angle-of-attack ``α``, side-slip angle ``β``, and a quasi-steady rotation vector ``\\Omega``.
"""
Freestream(V, α_deg, β_deg, Ω :: AbstractVector{T}) where T <: Real = Freestream{T}(V, α_deg, β_deg, Ω)

# TODO: ForwardDiff testing
# Freestream(V :: Real, α_deg, β_deg, Ω :: AbstractVector{<: Real}) = Freestream{Float64}(V, α_deg, β_deg, Ω)

"""
    Freestream(U, Ω)

A Freestream flow in Cartesian coordinates with vector ``U`` and quasi-steady rotation vector ``\\Omega``.
"""
Freestream(U :: AbstractVector{T}, Ω :: AbstractVector{T}) where T <: Real =
    let (V, α, β) = cartesian_to_freestream(U); Freestream{T}(V, α, β, Ω) end

"""
    freestream_to_cartesian(r, θ, φ)

Convert freestream flow (spherical polar) coordinates to Cartesian coordinates.
"""
freestream_to_cartesian(r, θ, φ) = r .* SVector(cos(θ) * cos(φ), -sin(φ), sin(θ) * cos(φ))

"""
    cartesian_to_freestream(U)

Convert Cartesian coordinates to freestream (spherical polar) flow coordinates.
"""
cartesian_to_freestream(U) = SVector(norm(U), -atand(U[3], U[1]), -atand(U[2], √(U[1]^2 + U[3]^2)))

"""
    velocity(freestream :: Freestream)

Compute the velocity of a `Freestream`.
"""
velocity(freestream :: Freestream) = freestream_to_cartesian(freestream.V, freestream.alpha, freestream.beta)

"""
    aircraft_velocity(freestream :: Freestream)

Compute the velocity of Freestream in the aircraft reference frame.
"""
aircraft_velocity(freestream :: Freestream) = -velocity(freestream)

end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
module QOptCtrl

using LinearAlgebra, Plots, Optim

include("basic_funcs.jl")
include("cost_funcs.jl")
include("optimization.jl")
include("propagator.jl")
include("visualization.jl")

export pauli_I, pauli_X, pauli_Y, pauli_Z, ket_0, ket_1, commutator, expval, outer_prod, tr_norm
export performance_func, C1_cost, C2_cost, density_op_fid, gate_infid, state_infid
export Basic_GRAPE, GRAPE, C1_grad, OPTIMIZE
export evolve, density_op, back_target_op, evolve_back
export visualise_pulses

end # module


text:
"""
    RootTaskREPL creates a Julia REPL where the REPL backend runs on the
    current Task and the frontend runs on a new task. This allows the
    backend to execute on the root Task if used at startup.

    If initialized on the root Task, the module will automatically start
    the REPL.

    This is useful for packages that must be initialized on the root Task
    on some platforms (e.g. JavaCall)

    Usage:
    julia -e "import RootTaskREPL"

    julia-root> Base.roottask === Base.current_task()
    true

    Advanced Usage:
    ROOTTASKREPL_ON_INIT=0 julia -e \\
    "import RootTaskREPL; RootTaskREPL.run_main_repl()"
"""
module RootTaskREPL
    include("repl.jl")

    ROOTTASKREPL_ON_INIT = true

    function __init__()
        global ROOTTASKREPL_ON_INIT = get(ENV,"ROOTTASKREPL_ON_INIT","1") in ("1","yes")
        if ROOTTASKREPL_ON_INIT && Base.roottask === Base.current_task()
           run_main_repl()
        end
    end

end # module


text:
function cable_transform(y, z)
    v1 = [0.0, 0.0, 1.0]
    v2 = y[1:3,1] - z[1:3,1]
    normalize!(v2)
    ax = cross(v1, v2)
    ang = acos(v1'*v2)
    R = AngleAxis(ang, ax...)

    if any(isnan.(R))
        R = I
    else
        nothing
    end

    compose(Translation(z), LinearMap(R))
end

function default_background!(vis)
    setvisible!(vis["/Background"], true)
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setvisible!(vis["/Axes"], false)
end


text:
using ReplMaker, Markdown

let
    keys() = println(join(replset.keys, ", "))

    global dashboard = false

    global replset = LogicalCombo()

    global setreplset!(x) = replset = x

    preserver = missing

    global activehistory = History()
    global setactivehistory!(x) = activehistory = x
    global pushactivehistory!(x) = push!(activehistory, x)

    global sessionhistory = History()
    global setsessionhistory!(x) = sessionhistory = x
    global pushsessionhistory!(x) = push!(sessionhistory, x)
    global popsessionhistory!() = pop!(sessionhistory)

    verboseall = true

    global replerror = false
    global replthrow(x) = (println(x) ; replerror = true)
    global returnreplerror() = replerror

    global preserver = missing
    global setpreserver!(x) = preserver = x

    global function abstractlogic(replinput; returnactive = false, verbose = true)
        replerror = false
        userinput = replinput |> strip |> tounicode

        userinput =  replace(userinput, r"^abstractlogic>"=>"") |> strip |> string

        verbose = verboseall & verbose & !occursin("[silent]", userinput)
        userinput = replace(userinput, "[silent]"=>"") |> strip |> string

        occursin("[clear]", userinput) && ALclear(verbose = verbose)
        userinput = replace(userinput, "[clear]"=>"") |> strip |> string

        if occursin(";;", userinput)
          for v in split(userinput, ";;");
              abstractlogic(v, verbose=verbose)
          end
          # userinput = "" # replace with return?
          returnactive && return replset
          return
        end

        # println("User input {$userinput}")

        if userinput == ""                                nothing(verbose = verbose)
        elseif (userinput[1] == '#')                      nothing(verbose = false)
        elseif occursin(r"⊥|dependenton|independentof", userinput)
            dependenton(userinput, replset, verbose = verbose)
        elseif occursin(r"[⊂⊃⊅⊄⋂⋔]|subset|superset|intersect|disjoint", userinput)
            setcompare(userinput, replset, verbose = verbose)
        elseif occursin(r"^t(est)?\(.*\)", userinput)                  testcall(userinput, verbose = verbose)
        elseif occursin(r"^(\?|help)", userinput)         help(userinput)
        elseif occursin(r"^show$", userinput)             ALshow(verbose = verbose)
        elseif occursin(r"^showall$", userinput)          ALshow(verbose = verbose, n = nfeasible(replset))
        elseif userinput ∈ ["back", "b"]                  ALback()
        elseif occursin(r"compare ", userinput)           ALcompare(userinput)
        elseif userinput ∈ ["discover", "d"]              discover(LogicalCombo())
        elseif occursin(r"^export( as){0,1} ", userinput) ALexport(userinput)
        elseif userinput ∈ ["next", "n", "f"]             ALnext(verbose = verbose)
        elseif occursin(r"^import ", userinput)           ALimport(userinput)
        elseif occursin(r"^dash(board)?$", userinput)     dashboard = !dashboard
        elseif userinput ∈ ["history", "h"]               Alhistory()
        elseif userinput ∈ ["History", "H"]               Alhistory(sessionprint=true)
        elseif userinput ∈ ["clearall", "Clear"]          ALClear(verbose = verbose)
        elseif occursin(r"^clear[ ]*$", userinput)        ALclear(verbose = verbose)
        elseif occursin(r"^range", userinput)             ALrange(userinput)
        elseif userinput ∈ ["keys", "k"]                  keys()
        elseif userinput == "preserve"                    ALpreserve()
        elseif userinput == "restore"                     ALrestore()
        elseif userinput == "silence"                     verboseall = false
        elseif userinput == "noisy"                       verboseall = true

        elseif occursin(r"^(prove|all|check|any|✓)", userinput) ALcheck(userinput)
        elseif occursin(r"^search", userinput)            ALsearch(userinput)
        else                                              ALparse(userinput, replset, verbose = verbose)
        end
        returnactive && return replset
        nothing
    end

end
"""
    abstractlogic(; returnactive = false)

Call the REPL from Julia. Setting returnactive to *true* returns the last active
logicset.

```julia
julia> abstractlogic("a, b, c in 1:3")
a, b, c ∈ 1:3            feasible outcomes 27 ✓          :1 1 3

julia> abstractlogic("a == b")
a == b                   feasible outcomes 9 ✓           :3 3 3

julia> abstractlogic("a > c")
a > c                    feasible outcomes 3 ✓           :3 3 1

julia> abstractlogic("c != 1")
c != 1                   feasible outcomes 1 ✓✓          :3 3 2

```
"""
abstractlogic = abstractlogic

function parse_to_expr(s)
   abstractlogic(s)

   # if dashboard
   #   # println("\nCommand Lists:" * string(replset.commands[end])
   #   # for v in showcommandlist()
   #   #     println(v)
   #   # end
   #   #
   #   # println("\nCommand Location:" * string(showcmdlocation()))
   #   # for v in showcommandhistory()
   #   #     println(v)
   #   # end
   #   #
   #   # println("\nLogicSet Lists:" * string(showsetlocation()))
   #   # for v in returnlogicset()
   #   #     println(join(collect(string(v))[1:(min(200,end))]))
   #   # end
   #   #
   # end
   nothing
end

"""
    replerror

Reports the error value of the most recently run repl command. Is defined as a
global in the `let block`.
"""
replerror = replerror

"""
    replset

Reports the value of the most recently run repl command. Is defined as a
global in the `let block`.

#### Example
```
abstractlogic> clear;; a, b in 1:2; a>b
Activeset Already Empty

a, b ∈ 1:2               Feasible Outcomes: 4    Perceived Outcomes: 4 ✓         :2 2
a>b                      Feasible Outcomes: 1    Perceived Outcomes: 1 ✓✓        :2 1

julia> replset |> nfeasible
1
```
"""
replset = replset


text:
# Define a minimal testing system with better reporting than the standard Julia one...
# This is taken from GodelTest; should be its own autotest package??

using Base.Test

module Minitest

type TestReporter
  start_time::Float64
  num_successes::Integer
  failures
  errors

  function TestReporter()
    new(time(), 0, Any[], Any[])
  end
end

set_start_time(tr::TestReporter) = tr.start_time = time()

num_failures(tr::TestReporter) = length(tr.failures)
num_errors(tr::TestReporter) = length(tr.errors)
num_assertions(tr::TestReporter) = tr.num_successes + num_failures(tr) + num_errors(tr)

pluralize(str::AbstractString, num) = (num == 1) ? str : (str * "s")

function report(tr::TestReporter)
  elapsed = time() - tr.start_time

  for i in 1:length(tr.failures)
    println(i, ") Failure:\n", tr.failures[i].expr)
  end

  for i in 1:length(tr.errors)
    println(i, ") Error:\n", tr.errors[i].expr)
    rethrow(tr.errors[i])
  end

  numassertions = num_assertions(tr)
  assertions_per_sec = numassertions/elapsed
  println("\n\n",
    @sprintf("Finished in %.2f seconds (%.1f assertions/sec)", elapsed, assertions_per_sec)
  )

  assertions_string = numassertions == 1 ? "assertion" : "assertions"

  numfailures = num_failures(tr)
  numerrors = num_errors(tr)

  println("\n",
    numassertions, pluralize(" assertion", numassertions), ", ",
    numfailures, pluralize(" failure", numfailures), ", ",
    numerrors, pluralize(" error", numerrors)
  )

  if (numfailures + numerrors) > 0
    exit(-1)
  else
    exit(0)
  end
end

TR__ = TestReporter()

CountsPerLine = 50
global Minitest_counts = 0

function log_result(indicator)
  global Minitest_counts
  print(indicator)
  Minitest_counts += 1
  if mod(Minitest_counts, CountsPerLine) == 0
    print(" $(Minitest_counts)\n")
  end
  flush(STDOUT)
end

# Include the file but also print to stdout and indent with spaces so that
# the log looks nice.
function include(filename)
  println("\n\nLoading $(filename):")
  global Minitest_counts
  spaces = map((x) -> " ", 1:mod(Minitest_counts, CountsPerLine))
  print(join(spaces))
  Main.include(filename)
end

function minitest_handler(r::Test.Success)
  TR__.num_successes += 1
  log_result(".")
end

function minitest_handler(r::Test.Failure)
  push!(TR__.failures, r)
  log_result("F")
end

function minitest_handler(r::Test.Error)
  push!(TR__.errors, r)
  log_result("E")
  #rethrow(r)
end

function do_tests(block)
  set_start_time(TR__)
  Base.Test.with_handler(block, minitest_handler)
  report(TR__)
end

end

macro repeatedly(expr)
  num_reps = isdefined(:NumTestRepetitions) ? NumTestRepetitions : 10
  quote
    for i in 1:$(num_reps)
      $(expr)
    end
  end
end


text:
x :: Type[Tag()]
x = Type

function f()

end


text:
"""
    find_word(
        graph::Dict{Int64,Array{String,1}},
        visited::Array{Bool,2},
        i::Int64,
        j::Int64,
        size_v::Int64,
        size_h::Int64,
        current_word::String,
        reference_words::Array{String,1},
        result::Array{String,1},
    )

Based on the Depth First Traversal algorithm, words will be find for a current char by
going through the graph up and down and keep track of the visited nodes. Important is that
the travelling happens in both direction up and down and left and right. Everythign will be
stacked in the visited list. If no word will be found the `current_word` has to be deleted.

# Arguments
- `graph::Dict{Int64,Array{String,1}}`: Graph of the connected nodes of chars, which can build the words
- `visited::Array{Bool,2}`: List of the visited chars in the graph
- `i::Int64`: current row
- `j::Int64`: current col
- `size_v::Int64`: total size of rows
- `size_h::Int64`: total size of cols
- `current_word::String`: current joint word of chars
- `reference_words::Array{String,1}`: Reference words to search for
- `result::Array{String,1}`: List of the found words in the graph
"""
function find_word(
    graph::Dict{Int64,Array{String,1}},
    visited::Array{Bool,2},
    i::Int64,
    j::Int64,
    size_v::Int64,
    size_h::Int64,
    current_word::String,
    reference_words::Array{String,1},
    result::Array{String,1},
)
    # Mark current vertex (CHAR) as visited
    visited[i, j] = true

    # Add visit vertex to current_word, which will be set to empty string by every new call
    current_word *= graph[i][j]

    #  If current_word is present in reference_words, then return it and stop
    if current_word in reference_words && !(current_word in result)
        push!(result, String(current_word))
    else
        # Traverse adjacent cells of graph
        for k in 1:i+1
            if k <= size_v
                for l in 1:j+1
                    if l <= size_h
                        if !visited[k, l]
                            find_word(
                                graph,
                                visited,
                                k,
                                l,
                                size_v,
                                size_h,
                                current_word,
                                reference_words,
                                result,
                            )
                        end
                    end
                end
            end
        end

        # Erase current character from string
        current_word = chop(current_word, tail = length(current_word) - 1)
        # Mark current vertex (CHAR) as not visited anymore
        visited[i, j] = false

    end
    return result

end

"""
    boogle_word_check(
        graph::Dict{Int64,Array{String,1}},
        reference_words::Array{String,1},
    )

For finding words (`reference_words`) in a field of chars, the boogle word check algorithm
goes for every single char up and down to see if the sum of the chars build a word
contained in th refernce word list. For this porpose the  Depth First Traversal algorithm
in function find_word is used.

# Arguments
- `graph::Dict{Int64,Array{String,1}}`: Graph of the connected nodes of chars, which can build the words
- `reference_words::Array{String,1}`: Reference words to search for

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> word_list = ["GEEKS", "FOR", "QUIZ", "GO"]
julia> graph_boogle = Dict(
        1 => ["G", "I", "Z"],
        2 => ["U", "E", "K"],
        3 => ["Q", "S", "E"],
        4 => ["D", "O", "P"],
        5 => ["F", "O", "R"],
        )
julia> ClassicAlgorithmsCollections.boogle_word_check(graph_boogle, word_list)
["GEEKS", "QUIZ", "FOR"]
```
"""
function boogle_word_check(
    graph::Dict{Int64,Array{String,1}},
    reference_words::Array{String,1},
)

    size_v = length(keys(graph))
    size_h = length(collect(values(graph))[1])
    #Create an array of visited nodes
    visited = zeros(Bool, (size_v, size_h))

    # Create a array list for the results
    result = String[]

    init_str = string()
    for i in 1:size_v
        for j in 1:size_h
            result = find_word(
                graph,
                visited,
                i,
                j,
                size_v,
                size_h,
                init_str,
                reference_words,
                result,
            )

        end
    end
    return result
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# adapted from http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw2a.py

import GDALUtils; const GU = GDALUtils

GU.registerdrivers() do
    open("ospy/data2/ut_counties.txt", "r") do file
        GU.create("tmp/hw2a.shp", "ESRI Shapefile") do output
            # create the layer for the output data source
            layer = GU.createlayer(output, "hw2a", GDAL.wkbPolygon)

            # create a field for the county name
            fielddefn = GU.createfielddefn("name", GDAL.OFTString)
            GU.setwidth(fielddefn, 30)

            # add the field to the shapefile
            GU.createfield(layer, fielddefn)

            # get the FeatureDefn for the shapefile
            featuredefn = GU.getlayerdefn(layer)

            # loop through the lines in the text file
            for line in readlines(file)
                # create an empty ring geometry
                ring = GU.creategeom(GDAL.wkbLinearRing)
                # split the line on colons to get county name and
                # a string with coordinates
                (name, coords) = split(line, ":")
                # split the coords on commas to get a list of x,y pairs
                coordlist = split(coords, ",")
                # loop through the list of coordinates
                for coord in coordlist
                    # split the x,y pair on spaces to get the individual x and y values
                    xy = split(coord)
                    # add the vertex to the ring
                    GU.addpoint(ring, parse(Float64, xy[1]),
                                      parse(Float64, xy[2]))
                end
                # now that we've looped through all of the coordinates, create a polygon
                poly = GU.creategeom(GDAL.wkbPolygon)
                GU.addgeom(poly, ring)

                # create a new feature and set its geometry and attribute
                GU.createfeature(featuredefn) do feature
                    GU.setgeom(feature, poly)
                    GU.setfield(feature, 0, name)

                    # add the feature to the shapefile
                    GU.createfeature(layer, feature)
                end
            end
        end
    end
end


text:
export predict

function predict(d::ConditionalModel, post::Vector{NamedTuple{N,T}}) where {N,T}
    args = argvals(d)
    m = d.model
    pred = predictive(m, keys(post[1])...)
    map(nt -> rand(pred(merge(args,nt))), post)
end

function predict(m::Model, post::Vector{NamedTuple{N,T}}) where {N,T}
    pred = predictive(m, keys(post[1])...)
    map(nt -> rand(pred(nt)), post)
end

# TODO: These don't yet work properly t on particles

function predict(d::ConditionalModel, post::NamedTuple{N,T}) where {N,T}
    args = argvals(d)
    m = Model(d)
    pred = predictive(m, keys(post)...)
    rand(pred(merge(args,post)))
end

function predict(m::Model, post::NamedTuple{N,T}) where {N,T}
    pred = predictive(m, keys(post)...)
    rand(pred(post))
end

predict(m::Model; kwargs...) = predict(m,(;kwargs...))

predict(d,x) = x

function predict(d, s::AbstractVector)
    [predict(d, sj) for sj in s]
end


text:
@doc """
      Euler Method
      Supported by:
      -PhD Kelvyn B. Sánchez kelvyn.baruc@gmail.com
 """ ->
function euler(f, tspan, y0; h = 1)
t = tspan[1]:h:tspan[end]
n = length(t)

      if t[n] < tspan[end]
            t = [t, tspan[end]]
            n += 1
      end

y1 = f(t[1],y0)

eq = length(y1) # Test the number of equations in f
y = zeros(Int(n),eq)
y[1,:] = y0

for i=1:n-1 # Solve for each x
      for j=1:eq # Solve for each equation in f
            y1 = f(t[i],y[i,:])
            y[i+1,j] = y[i,j] +  y1[j] * (t[i+1] - t[i])

      end
end
return t, y
end


text:
using Random

"""
Example Use:
> genObject("object1", env)
"""

env = Dict(["custom_types" => Dict([
            "Object1" => [],
            "Object2" => [],
            "Object3" => [],
            ]),
            "variables" => Dict([
              "object1" => "Object_Object2",
              "object2" => "Object_Object3",
              "objectlist1" => "ObjectList_Object3",
            ])])

function genUpdateRule(var, environment; p=0.7)
  if environment["variables"][var] == "Int"
    genInt(environment)
  elseif environment["variables"][var] == "Bool"
    genBoolUpdateRule(var, environment)
  elseif occursin("Object_", environment["variables"][var])
    genObjectUpdateRule(var, environment, p=p)
  else
    genObjectListUpdateRule(var, environment, p=p)
  end
end

# -----begin object generator + helper functions ----- #
function genObject(environment; p=0.9)
  object = genObjectName(environment)
  genObjectUpdateRule(object, environment, p=p)
end

function genObjectName(environment)
  objects = filter(var -> occursin("Object_", environment["variables"][var]), collect(keys(environment["variables"])))
  objects[rand(1:length(objects))]
end

function genObjectUpdateRule(object, environment; p=0.0)
  prob = rand()
  if prob < p
    if object == "obj"
      "$(object)"
    else
      "(prev $(object))"
    end
  else

    choices = [
      # "(moveLeft $(object))",
      # ("moveLeft", [:(genObjectUpdateRule($(object), $(environment), p=0.9))]),
      # ("moveRight", [:(genObjectUpdateRule($(object), $(environment), p=0.9))]),
      # ("moveUp", [:(genObjectUpdateRule($(object), $(environment), p=0.9))]),
      # ("moveDown", [:(genObjectUpdateRule($(object), $(environment), p=0.9))]),
      # ("moveNoCollision", [:(genObjectUpdateRule($(object), $(environment))), :(genPosition($(environment)))]),
      # ("nextLiquid", [:(genObjectUpdateRule($(object), $(environment)))]),
      # ("nextSolid", [:(genObjectUpdateRule($(object), $(environment)))]),
      # ("rotate", [:(genObjectUpdateRule($(object), $(environment)))]),
      # ("rotateNoCollision", [:(genObjectUpdateRule($(object), $(environment)))]),
    ]
    other_types = filter(type -> type != environment["variables"][object], collect(keys(environment["custom_types"])))
    if length(other_types) > 0
      push!(choices, "(move $(object) (unitVector $(object) (closest $(object) $(split(rand(other_types), "_")[2]))))")
    end
    if choices == []
      choice = object
    else
      choice = choices[rand(1:length(choices))]
    end
    # "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
    choice
  end
end

# ----- end object generator + helper functions ----- #

# ----- Int generator ----- #
function genInt(environment)
  int_vars = map(v -> "(prev $(v))", filter(var -> environment["variables"][var] == "Int", collect(keys(environment["variables"]))))
  choices = [ #= fieldsFromCustomTypes("Int", environment)..., =# collect(1:5)..., int_vars...]
  choice = rand(choices)
  if (choice isa String) || (choice isa Int)
    choice
  else
    @show "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
    "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
  end
end

# ----- Bool generator ----- #
function genBool(environment)
  choices = [
    ("clicked", []),
    ("clicked", [:(genPosition($(environment)))]),
    ("left", []),
    ("right", []),
    ("up", []),
    ("down", []),
    ("true", []), # TODO: add not, or, and -- need to be able to specify prior probabilities
  ]
  if length(filter(var -> occursin("Object_", environment["variables"][var]), collect(keys(environment["variables"])))) > 0
    push!(choices, [("clicked", [:(genObject($(environment), p=1.0))]),
                    ("intersects", [:(genObject($(environment), p=1.0)), :(genObject($(environment), p=1.0))]),
                    ("isWithinBounds", [:(genObject($(environment)))])
                   ]...)
  end

  bool_vars = map(v -> "(prev $(v))", filter(var -> environment["variables"][var] == "Bool", collect(keys(environment["variables"]))))
  foreach(var -> push!(choices, (var, [])), bool_vars)

  push!(choices, fieldsFromCustomTypes("Bool", environment)...)

  choice = choices[rand(1:length(choices))]
  if (length(choice[2]) == 0)
    "$(choice[1])"
  else
    "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
  end

end

function genBoolUpdateRule(bool, environment)
  bool_val = "(prev $(bool))"
  rand(["(! $(bool_val))", bool_val])
end

# ----- Position generator ----- #
function genPosition(environment)
  choices = [
    ("Position", [:(genInt($(environment))), :(genInt($(environment)))]),
    ("displacement", [:(genPosition($(environment))), :(genPosition($(environment)))]),
    ("unitVector", [:(genPosition($(environment)))]),
    ("uniformChoice", [:(genPositionList($(environment)))])
  ]
  # if object constants exist, add support for (.. obj origin)
  if length(filter(var -> occursin("Object_", environment["variables"][var]), collect(keys(environment["variables"])))) > 0
    push!(choices, ("..", [:(genObject($(environment), p=1.0)), :(String("origin"))]))
  end

  choice = choices[rand(1:length(choices))]
  "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
end

# ----- Click generator ----- #
function genClick(environment)
  options = [
    "click"
  ]
  choice = choices[rand(1:length(choices))]
  "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
end

# ----- Position List generator ----- #
function genPositionList(environment)
  choices = [
    ("randomPositions", ["GRID_SIZE", :(genInt($(environment)))])
  ]
  if length(filter(var -> occursin("Object_", environment["variables"][var]), collect(keys(environment["variables"])))) > 0
    push!(choices, ("adjPositions", [:(genObject($(environment)))]))
  end
  choice = choices[rand(1:length(choices))]
  "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
end

# ----- begin object list generator + helper functions ----- #
function genObjectList(environment)
  object_list = genObjectListName(environment)
  genObjectListUpdateRule(object_list, environment)
end

function genObjectListName(environment)
  object_lists = filter(var -> occursin("ObjectList_", environment["variables"][var]), collect(keys(environment["variables"])))
  object_lists[rand(1:length(object_lists))]
end

function genObjectConstructor(type, environment)
  new_type = occursin("_", type) ? type : string("Object_", type)
  constructor = map(tuple -> Meta.parse("gen$(tuple[2])($(environment))"), environment["custom_types"][new_type])
  push!(constructor, :(genPosition($(environment))))
  "($(type) $(join(map(eval, constructor), " ")))"
end

function genObject(type, environment, p=0.9)
  objects_with_type = filter(var -> environment["variables"][var] == type, collect(keys(environment["variables"])))
  prob = rand()
  if (prob < p) && length(objects_with_type) != 0
    rand(objects_with_type)
  else
    constructor = genObjectConstructor(type, environment)
    constructor
  end
end

function genObjectListUpdateRule(object_list, environment; p=0.7)
  prob = rand()
  if prob < p
    "(prev $(object_list))"
  else
    choices = [
      ("addObj",
        [:(genObjectListUpdateRule($(object_list), $(environment))),
         :(genObjectConstructor($(String(split(environment["variables"][object_list], "_")[end])), $(environment))),
        ]
      ),
      ("updateObj",
        [:(genObjectListUpdateRule($(object_list), $(environment))),
         :(genLambda($(environment)))
        ]),
    ]
    choice = choices[rand(1:length(choices))]
    "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
  end
end

function genLambda(environment)
  choice = ("-->", [:(String("obj")), :(genObjectUpdateRule("obj", $(environment)))])
  "($(choice[1]) $(join(map(eval, choice[2]), " ")))"
end
# ----- end object list generator + helper functions ----- #

# ----- string generator ----- #
function genString(environment)
  colors = ["red", "yellow", "green", "blue"]
  color = colors[rand(1:length(colors))]
  """ "$(color)" """
end

# ----- helper functions ----- #

function fieldsFromCustomTypes(fieldtype::String, environment)
  branches = []
  types_with_field = filter(type -> fieldtype in map(tuple -> tuple[2], environment["custom_types"][type]), collect(keys(environment["custom_types"])))
  for type in types_with_field
    fieldnames = map(t -> t[1], filter(tuple -> tuple[2] == fieldtype, environment["custom_types"][type]))
    if length(filter(var -> environment["variables"][var] == type, collect(keys(environment["variables"])))) > 0
      foreach(fieldname -> push!(branches, ("..", [Meta.parse("genObject(\"$(split(type, "_")[end])\", $(environment))"), :(String($(fieldname)))])), fieldnames)
    end
  end
  branches
end


text:
function Grid=TriGrid(nx,ny,lx,ly,Case)
switch Case
  case 'Tri2'
    Grid=TriGrid2(nx,ny,lx,ly);
  case 'Tri4'
    Grid=TriGrid4(nx,ny,lx,ly);
end
Grid=RenumberingTri(Grid);
Grid=OrientationFaceTri(Grid);
Grid=OrdEdgeFacesTri(Grid);
Grid.Dim=3;
Grid.Type='Triangle';
Grid.nBar=[ 0   1   1
           -1   1   0];
end
function Grid=TriGrid2(nx,ny,lx,ly)
NumNodes=(nx+1)*(ny+1);
Nodes(1:NumNodes)=Node([0 0 0],0);
NodeNumber=1;
dx=lx/nx;
dy=ly/ny;
y=0;
for iy=1:ny+1
  x=0;
  for ix=1:nx+1
    Nodes(NodeNumber)=Node([x;y;0],NodeNumber);
    NodeNumber=NodeNumber+1;
    x=x+dx;
  end
  y=y+dy;
end
Grid.Nodes=Nodes;
Grid.NumNodes=NumNodes;

NumEdges=(nx+1)*ny+nx*(ny+1)+nx*ny;
Edges(1:NumEdges)=Edge([1 2],Grid,0);

N1=1;
N2=1+nx+1;
EdgeNumber=1;
for iy=1:ny
  for ix=1:nx+1
    Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber);
    EdgeNumber=EdgeNumber+1;
    N1=N1+1;
    N2=N2+1;
  end
end
N1=1;
N2=2;
for iy=1:ny+1
  for ix=1:nx
    Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber);
    EdgeNumber=EdgeNumber+1;
    N1=N1+1;
    N2=N2+1;
  end
  N1=N1+1;
  N2=N2+1;
end
N1=2;
N2=nx+2;
for iy=1:ny
  for ix=1:nx
    Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber);
    EdgeNumber=EdgeNumber+1;
    N1=N1+1;
    N2=N2+1;
  end
  N1=N1+1;
  N2=N2+1;
end
Grid.Edges=Edges;
Grid.NumEdges=NumEdges;

NumFaces=2*nx*ny;
Faces(1:NumFaces)=Face([0 0 0 0],Grid,0,'x');
E1=(nx+1)*ny+1;
E2=2;
E3=(nx+1)*ny+1+nx;
E4=1;
E5=(nx+1)*ny+nx*(ny+1)+1;
FaceNumber=1;
Type='o';
for iy=1:ny
  for ix=1:nx
    [Faces(FaceNumber),Grid]=Face([E1 E5 E4],Grid,FaceNumber,Type);
    FaceNumber=FaceNumber+1;
    [Faces(FaceNumber),Grid]=Face([E2 E3 E5],Grid,FaceNumber,Type);
    FaceNumber=FaceNumber+1;
    E1=E1+1;
    E2=E2+1;
    E3=E3+1;
    E4=E4+1;
    E5=E5+1;
  end
  E2=E2+1;
  E4=E4+1;
end
Grid.Faces=Faces;
Grid.NumFaces=NumFaces;
end

function Grid=TriGrid4(nx,ny,lx,ly)

  NumNodes=(nx+1)*(ny+1)+nx*ny;
  Nodes(1:NumNodes)=Node([0 0 0],0);
  NodeNumber=1;
  dx=lx/nx;
  dy=ly/ny;
  y=0;
  for iy=1:ny+1
    x=0;
    for ix=1:nx+1
      Nodes(NodeNumber)=Node([x;y;0],NodeNumber);
      NodeNumber=NodeNumber+1;
      x=x+dx;
    end
    y=y+dy;
  end
  y=0.5*dy;
  for iy=1:ny
    x=0.5*dx;
    for ix=1:nx
      Nodes(NodeNumber)=Node([x;y;0],NodeNumber);
      NodeNumber=NodeNumber+1;
      x=x+dx;
    end
    y=y+dy;
  end
  Grid.Nodes=Nodes;
  Grid.NumNodes=NumNodes;

  NumEdges=(nx+1)*ny+nx*(ny+1)+4*nx*ny;
  Edges(1:NumEdges)=Edge([1 2],Grid,0);

  N1=1;
  N2=1+nx+1;
  EdgeNumber=1;
  for iy=1:ny
    for ix=1:nx+1
      Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber);
      EdgeNumber=EdgeNumber+1;
      N1=N1+1;
      N2=N2+1;
    end
  end
  N1=1;
  N2=2;
  for iy=1:ny+1
    for ix=1:nx
      Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber);
      EdgeNumber=EdgeNumber+1;
      N1=N1+1;
      N2=N2+1;
    end
    N1=N1+1;
    N2=N2+1;
  end

  N1=1;
  N2=2;
  N3=nx+3;
  N4=nx+2;
  NC=(nx+1)*(ny+1)+1;
  for iy=1:ny
    for ix=1:nx
      Edges(EdgeNumber)=Edge([N1 NC],Grid,EdgeNumber);
      EdgeNumber=EdgeNumber+1;
      Edges(EdgeNumber)=Edge([N2 NC],Grid,EdgeNumber);
      EdgeNumber=EdgeNumber+1;
      Edges(EdgeNumber)=Edge([N3 NC],Grid,EdgeNumber);
      EdgeNumber=EdgeNumber+1;
      Edges(EdgeNumber)=Edge([N4 NC],Grid,EdgeNumber);
      EdgeNumber=EdgeNumber+1;
      N1=N1+1;
      N2=N2+1;
      N3=N3+1;
      N4=N4+1;
      NC=NC+1;
    end
    N1=N1+1;
    N2=N2+1;
    N3=N3+1;
    N4=N4+1;
  end
  Grid.Edges=Edges;
  Grid.NumEdges=NumEdges;

  NumFaces=4*nx*ny;
  Faces(1:NumFaces)=Face([0 0 0 0],Grid,0,'x');
  E1=(nx+1)*ny+1;
  E2=2;
  E3=(nx+1)*ny+1+nx;
  E4=1;
  EC1=(nx+1)*ny+nx*(ny+1)+1;
  EC2=(nx+1)*ny+nx*(ny+1)+2;
  EC3=(nx+1)*ny+nx*(ny+1)+3;
  EC4=(nx+1)*ny+nx*(ny+1)+4;
  FaceNumber=1;
  Type='o';
  for iy=1:ny
    for ix=1:nx
      [Faces(FaceNumber),Grid]=Face([EC1 E1 EC2],Grid,FaceNumber,Type);
      FaceNumber=FaceNumber+1;
      [Faces(FaceNumber),Grid]=Face([EC2 E2 EC3],Grid,FaceNumber,Type);
      FaceNumber=FaceNumber+1;
      [Faces(FaceNumber),Grid]=Face([EC3 E3 EC4],Grid,FaceNumber,Type);
      FaceNumber=FaceNumber+1;
      [Faces(FaceNumber),Grid]=Face([EC4 E4 EC1],Grid,FaceNumber,Type);
      FaceNumber=FaceNumber+1;
      E1=E1+1;
      E2=E2+1;
      E3=E3+1;
      E4=E4+1;
      EC1=EC1+4;
      EC2=EC2+4;
      EC3=EC3+4;
      EC4=EC4+4;
    end
    E2=E2+1;
    E4=E4+1;
  end
  Grid.Faces=Faces;
  Grid.NumFaces=NumFaces;

end

function [Grid] = RenumberingTri(Grid)
for iF=1:Grid.NumFaces
  Grid.Faces(iF)=RenumberingFace3(Grid.Faces(iF),Grid);
end
for iE=1:Grid.NumEdges
  Grid.Edges(iE)=PosEdgeInFace(Grid.Edges(iE),Grid);
end
end

function Edge=PosEdgeInFace(Edge,Grid)
Edge.FE=zeros(2,1);
for iF=1:size(Edge.F,2)
  F=Edge.F(iF);
  for iE=1:3
    if Edge.E==Grid.Faces(F).E(iE)
      Edge.FE(iF)=iE;
      break
    end
  end
end
end

function Face=RenumberingFace3(Face,Grid)
Face.N=sort(Face.N);
N=[Face.N Face.N(1)];
E=Face.E;
N2=N(1);
for iN=1:3
  N1=N2;
  N2=N(iN+1);
  for iE=1:3
    if (N1==Grid.Edges(E(iE)).N(1) && N2==Grid.Edges(E(iE)).N(2)) || ...
        (N1==Grid.Edges(E(iE)).N(2) && N2==Grid.Edges(E(iE)).N(1))
      Face.E(iN)=E(iE);
      break
    end
  end
end
end

function Grid=OrientationFaceTri(Grid)
for iF=1:Grid.NumFaces
  Face=Grid.Faces(iF);
  r=cross(Grid.Nodes(Face.N(1)).P-Grid.Nodes(Face.N(2)).P,Grid.Nodes(Face.N(3)).P-Grid.Nodes(Face.N(2)).P);
  if r(3)<0
    Grid.Faces(iF).Ori=1;
  else
    Grid.Faces(iF).Ori=-1;
  end
end
end

function Grid=OrdEdgeFacesTri(Grid)
for iE=1:Grid.NumEdges
  Edge=Grid.Edges(iE);
  if size(Edge.F,2)>1
    if iE==Grid.Faces(Edge.F(1)).E(3) && Grid.Faces(Edge.F(1)).Ori==1
      Grid.Edges(iE).F(1)=Edge.F(2);
      Grid.Edges(iE).F(2)=Edge.F(1);
      Grid.Edges(iE).FE(1)=Edge.FE(2);
      Grid.Edges(iE).FE(2)=Edge.FE(1);
    end
  end
end
end


text:
type CTreeChart <: JFinM_Charts.AbstractChart
    chartType::String
    cmd::String
    shape::String
    arrowhead::String
    emph1::Array{Int, 1}
    emphFillColor1::String
    emphShape1::String
    emph2::Array{Int, 1}
    emphFillColor2::String
    emphShape2::String
end

## define default values
CTreeChart() = CTreeChart("CTreeChart",
                          "dot",
                          "circle", "none",
                          [], "lawngreen", "circle",
                          [], "lawngreen", "circle")

##########
## GViz ##
##########

@doc doc"""
Transform data and chart into graphviz instance.
"""->
function GViz(tr::CTreeParRef, chrt::CTreeChart)
    dotCode = ""

    ## emph1 and emph2 sets need to be non-overlapping
    if !isempty(intersect(chrt.emph1, chrt.emph2))
        error("emphasized sets need to be non-overlapping")
    end

    nVars = size(tr.tree, 1)
    dotCode = string(dotCode, "digraph {\n")
    dotCode = string(dotCode, "node [shape=$(chrt.shape)];\n")
    dotCode = string(dotCode, "edge [arrowhead=$(chrt.arrowhead)];\n")
    for ii=1:nVars
        parNode = tr.tree[ii]
        if parNode > 0
            if ii in chrt.emph1
                dotCode = string(dotCode, "$ii [style=filled; fillcolor=$(chrt.emphFillColor1); shape=$(chrt.emphShape1)];")
            elseif ii in chrt.emph2
                dotCode = string(dotCode, "$ii [style=filled; fillcolor=$(chrt.emphFillColor2); shape=$(chrt.emphShape2)];")
            else

                end
            dotCode = string(dotCode, "$parNode -> $ii ;\n")
        end
    end
    dotCode = string(dotCode, "}")
    return GViz("CTreeChart", "dot", dotCode)
end

function GViz(tr::CTreePaths, chrt::CTreeChart)
    return GViz(convert(CTreeParRef, tr), chrt)
end


text:
using FFTViews
using Base.Test

function test_approx_eq_periodic(a::FFTView, b)
    for I in CartesianRange(indices(b))
        @test_approx_eq a[I-1] b[I]
    end
    nothing
end

function test_approx_eq_periodic(a::FFTView, b::FFTView)
    for I in CartesianRange(indices(b))
        @test_approx_eq a[I] b[I]
    end
    nothing
end

@testset "ambiguities" begin
    @test isempty(detect_ambiguities(FFTViews, Base, Core))
end

@testset "basics" begin
    a = FFTView{Float64,2}((5,7))
    @test indices(a) == (0:4, 0:6)
    @test eltype(a) == Float64
    a = FFTView{Float64}((5,7))
    @test indices(a) == (0:4, 0:6)
    @test eltype(a) == Float64
    @test_throws MethodError FFTView{Float64,3}((5,7))
    for i = 1:35
        a[i] = i
    end
    @test a[3,1] == 9
    @test a[:,0] == FFTView(collect(1:5))
    @test a[0,:] == FFTView(collect(1:5:35))
    @test a[1,0:7] == [2:5:35;2]
    @test a[2,[0,1,0,-1]] == [3,8,3,33]
    @test a[3,trues(9)] == [9,14,19,24,29,34,4,9,14]
    @test a[3,FFTView(trues(9))] == [4,9,14,19,24,29,34,4,9]
    b = similar(Array{Int}, indices(a))
    @test isa(b, FFTView)
    @test indices(b) == indices(a)
    @test eltype(b) == Int
    @test reshape(a, Val{2}) === a
    @test reshape(a, Val{1}) == FFTView(convert(Vector{Float64}, collect(1:35)))
    @test indices(reshape(a, Val{3})) == (0:4,0:6,0:0)
end

@testset "convolution-shift" begin
    for l in (8,9)
        a = zeros(l)
        v = FFTView(a)
        @test indices(v,1) == 0:l-1
        v[0] = 1
        p = rand(l)
        pfilt = ifft(fft(p).*fft(v))
        @test_approx_eq real(pfilt) p
        v[0] = 0
        v[-1] = 1
        pfilt = ifft(fft(p).*fft(v))
        @test_approx_eq real(pfilt) circshift(p, -1)
        v[-1] = 0
        v[+1] = 1
        pfilt = ifft(fft(p).*fft(v))
        @test_approx_eq real(pfilt) circshift(p, +1)
    end
    for l2 in (8,9), l1 in (8,9)
        a = zeros(l1,l2)
        v = FFTView(a)
        @test indices(v) == (0:l1-1, 0:l2-1)
        p = rand(l1,l2)
        for offset in ((0,0), (-1,0), (0,-1), (-1,-1),
                       (1,0), (0,1), (1,1), (1,-1), (-1,1),
                       (3,-5), (281,-14))
            fill!(a, 0)
            v[offset...] = 1
            pfilt = ifft(fft(p).*fft(v))
            @test_approx_eq real(pfilt) circshift(p, offset)
        end
    end
end

using OffsetArrays

@testset "convolution-offset" begin
    for l2 in (8,9), l1 in (8,9)
        a = OffsetArray(zeros(l1,l2), (-2,-3))
        v = FFTView(a)
        @test indices(v) == (-2:l1-3, -3:l2-4)
        p = rand(l1,l2)
        po = OffsetArray(copy(p), (5,-1))
        for offset in ((0,0), (-1,0), (0,-1), (-1,-1),
                       (1,0), (0,1), (1,1), (1,-1), (-1,1),
                       (3,-5), (281,-14))
            fill!(a, 0)
            v[offset...] = 1
            pfilt = ifft(fft(p).*fft(v))
            @test_approx_eq real(pfilt) circshift(p, offset)
            pofilt = ifft(fft(po).*fft(v))
            test_approx_eq_periodic(FFTView(real(pofilt)), circshift(po, offset))
            pfilt = irfft(rfft(p).*rfft(v), length(indices(v,1)))
            @test_approx_eq real(pfilt) circshift(p, offset)
            pofilt = irfft(rfft(po).*rfft(v), length(indices(v,1)))
            test_approx_eq_periodic(FFTView(real(pofilt)), circshift(po, offset))
            dims = (1,2)
            pfilt = ifft(fft(p, dims).*fft(v, dims), dims)
            @test_approx_eq real(pfilt) circshift(p, offset)
            pofilt = ifft(fft(po, dims).*fft(v, dims), dims)
            test_approx_eq_periodic(FFTView(real(pofilt)), circshift(po, offset))
            pfilt = irfft(rfft(p, dims).*rfft(v, dims), length(indices(v,1)), dims)
            @test_approx_eq real(pfilt) circshift(p, offset)
            pofilt = irfft(rfft(po, dims).*rfft(v, dims), length(indices(v,1)), dims)
            test_approx_eq_periodic(FFTView(real(pofilt)), circshift(po, offset))
        end
    end
end

@testset "vector indexing" begin
    v = FFTView(1:10)
    @test v[-10:15] == [1:10;1:10;1:6]
end

nothing


text:
using Distributions: Uniform, Normal
using StatsBase: mean, sample, standardize, ZScoreTransform
using Random
using LinearAlgebra
using MultivariateStats: ridge

function RocketFit(X::Array, y::AbstractArray, num_kernels::Int64=10000, normalise::Bool=true, random_state=nothing)
    n_timepoints, n_columns = size(X)  #n_timepoints should be max of sereis length in multidimentional series!
    _X = standardize(ZScoreTransform, X, dims=2)
    random_state = typeof(random_state) == Int64 ? random_state : nothing
    kernels = generate_kernels(n_timepoints, num_kernels, random_state)
    _X = apply_kernels(_X, kernels)
    Parameters = ridge(_X, y, 0)
    return Parameters, Kernels
end

function Predict_new(Xnew::Array, kernels, Parameters::Array)
    _Xnew = standardize(ZScoreTransform, Xnew, dims=2)
    _Xnew = apply_kernels(_Xnew, kernels)
    y_hat = _Xnew*Parameters[1:end-1] + Parameters[end]
    return y_hat
end

function generate_kernels(n_timepoints::Int64, num_kernels::Int64, seed::Union{Bool, Int64})
    if seed != false
        Random.seed!(seed)
    end
    candidate_lengths = Int32[7, 9, 11]
    lengths = sample(candidate_lengths, num_kernels)
    weights = zeros(Float64, sum(lengths))

    biases = zeros(Float64, num_kernels)
    dilations = zeros(Int32, num_kernels)
    paddings = zeros(Int32, num_kernels)
    a1 = 1
    for i=1:num_kernels
        _length = lengths[i]
        _weights = rand(Normal(),  _length)
        b1 = a1 + _length
        weights[a1:b1-1] = _weights .- mean(_weights)
        biases[i] = rand(Uniform(-1, 1))
        dilation = floor(Int32, 2^( rand( Uniform(0, log2((n_timepoints-1)/(_length-1)) ))))
        dilations[i] = dilation
        padding = rand((0,1)) == 1 ? ( ((_length - 1) * dilation) ÷ 2) : 0      # check //
        paddings[i] = padding
        a1 = b1
    end
    return weights, lengths, biases, dilations, paddings
end

function apply_kernel(X::AbstractArray, weights::AbstractArray, length_l::Int32,
                                         bias::Float64, dilation::Int32, padding::Int32) # change name if used as length above
    n_timepoints = length(X)
    output_length = (n_timepoints + (2 * padding)) - ((length_l - 1) * dilation)
    _ppv = 0
    _max = -maxintfloat(Float64)
    endl = (n_timepoints + padding) - ((length_l - 1) * dilation)
    for i = -padding:(endl-1)
        _sum = bias
        index = i + 1   # This is change, check for dilation > 1
        for j =1:(length_l)
            if index > 0 && index <= n_timepoints    # and in julia
                _sum = _sum + weights[j] * X[index]
            end
            index = index + dilation
        end
        if _sum > _max
            _max = _sum
        end
        if _sum > 0
            _ppv += 1
        end
    end
   return _ppv/output_length, _max
end

function apply_kernels(X::Array, kernels)
    weights, lengths, biases, dilations, paddings = kernels
    #n_columns = 1
    n_instances, _ = size(X)
    num_kernels = length(lengths)
    _X = zeros(Float64, n_instances, (num_kernels * 2))  # 2 features per kernel
    for i = 1:n_instances
        a1 = 1  # for weights
        a2 = 1  # for features
        for j =1:num_kernels
            b1 = a1 + lengths[j]
            b2 = a2 + 2
            _X[i, a2:b2-1] .= apply_kernel(X[i, :], weights[a1:b1-1],
                                      lengths[j], biases[j], dilations[j], paddings[j]) #check the =

            a1 = b1
            a2 = b2
        end
    end
    return _X
end

#=
These methods will be used for multivariate problems.
function _generate_kernels(n_timepoints, num_kernels, n_columns, seed)
    if seed != false
        Random.seed!(seed)
    end

    candidate_lengths = Int32[7, 9, 11]
    lengths = sample(candidate_lengths, num_kernels)
    num_channel_indices = zeros(Int32, num_kernels)

    @inbounds for i=1:num_kernels
        limit = min(n_columns, lengths[i])
        num_channel_indices[i] = floor(Int32, 2^(rand(Uniform(0, log2(limit + 1)))))
    end

    channel_indices = zeros(Int32, floor(Int32, sum(num_channel_indices)))

    weights = zeros(Float64, floor(Int32, dot(lengths, num_channel_indices)))

    biases = zeros(Float64, num_kernels)
    dilations = zeros(Int32, num_kernels)
    paddings = zeros(Int32, num_kernels)

    a1 = 1  # for weights
    a2 = 1  # for channel_indices

    @inbounds for i=1:num_kernels

        _length = lengths[i]
        _num_channel_indices = num_channel_indices[i]

        _weights = rand(Normal(), _num_channel_indices * _length)

        b1 = a1 + (_num_channel_indices * _length)
        b2 = a2 + _num_channel_indices

        a3 = 1  # Check index
        @inbounds for _ = 1:_num_channel_indices
            b3 = a3 + _length
            _weights[a3:b3-1] = _weights[a3:b3-1] .- mean(_weights[a3:b3-1])
            a3 = b3
        end
        weights[a1:b1-1] = _weights

        channel_indices[a2:b2-1] = sample([0:(n_columns-1);], _num_channel_indices, replace=false)

        biases[i] = rand(Uniform(-1, 1))

        dilation = floor(Int32, 2^( rand( Uniform(0, log2((n_timepoints-1)/(_length-1)) ))))

        dilations[i] = dilation

        padding = rand((0,1)) == 1 ? ( ((_length - 1) * dilation) ÷ 2) : 0      # check //
        paddings[i] = padding

        a1 = b1
        a2 = b2
    end

    return (
        weights,
        lengths,
        biases,
        dilations,
        paddings,
        num_channel_indices,
        channel_indices,
       )
end

function _apply_kernel_univariate(X, weights, length_l, bias, dilation, padding) # change name if used as length above
    n_timepoints = length(X)

    output_length = (n_timepoints + (2 * padding)) - ((length_l - 1) * dilation)

    _ppv = 0
    _max = -maxintfloat(Float64)

    endl = (n_timepoints + padding) - ((length_l - 1) * dilation)

    @inbounds for i = -padding:(endl-1)

        _sum = bias

        index = i

        @inbounds for j =1:(length_l)

            if index > -1 && index < n_timepoints    # and in julia
                _sum = _sum + weights[j] * X[index]
            end
            index = index + dilation
        end
        if _sum > _max
            _max = _sum
        end
        if _sum > 0
            _ppv += 1
        end
    end
    return _ppv/output_length, _max
end

function _apply_kernel_multivariate(X, weights, length_l, bias, dilation, padding, num_channel_indices, channel_indices)
    n_timepoints = size(X)

    output_length = (n_timepoints + (2 * padding)) - ((length_l - 1) * dilation)

    _ppv = 0
    _max = -maxintfloat(Float64)

    endl = (n_timepoints + padding) - ((length_l - 1) * dilation)

    @inbounds for i = -padding:(endl-1)

        _sum = bias

        index = i

        @inbounds for j =1:(length_l)

            if index > -1 && index < n_timepoints # and in julia

                @inbounds for k = 1:num_channel_indices
                    _sum = _sum + weights[k, j] * X[channel_indices[k], index]
                end
            index = index + dilation
            end
        end
        if _sum > _max
            _max = _sum
        end
        if _sum > 0
            _ppv += 1
        end
    end
    return _ppv / output_length, _max
end

function _apply_kernels(X, kernels)

    weights, lengths, biases, dilations, paddings, num_channel_indices, channel_indices = kernels
    n_columns = 1
    n_instances, _ = size(X)
    num_kernels = length(lengths)

    _X = zeros(Float64, n_instances, (num_kernels * 2))  # 2 features per kernel

    for i = 1:n_instances   # check waht is Prange!!

        a1 = 1  # for weights
        a2 = 1  # for channel_indices
        a3 = 1  # for features

        for j =1:num_kernels

            b1 = a1 + num_channel_indices[j] * lengths[j]
            b2 = a2 + num_channel_indices[j]
            b3 = a3 + 2

            if num_channel_indices[j] == 1

                _X[i, a3:b3-1] = _apply_kernel_univariate(
                    X[i, channel_indices[a2-1]],
                    weights[a1:b1-1],
                    lengths[j],
                    biases[j],
                    dilations[j],
                    paddings[j],
                )

            else

                _weights = reshape(weights[a1:b1-1], (num_channel_indices[j], lengths[j]))

                _X[i, a3:b3-1] = _apply_kernel_multivariate(
                    X[i, :],
                    _weights,
                    lengths[j],
                    biases[j],
                    dilations[j],
                    paddings[j],
                    num_channel_indices[j],
                    channel_indices[a2:b2-1]
                )
            end
            a1 = b1
            a2 = b2
            a3 = b3
        end
    end
    return _X
end
=#


text:
"""
Contains stubs for advective, diffusive, and nondiffusive courant number
calculations to be used in [`DGmethods.courant`](@ref). Models should provide
concrete implementations if they wish to use these functions.
"""
module Courant

export advective_courant, diffusive_courant, nondiffusive_courant

function advective_courant end

function nondiffusive_courant end

function diffusive_courant end

function viscous_courant end

end


text:
module Pardiso

if !isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
    error("""please run Pkg.build("Pardiso") before loading the package""")
end

using Compat

import Compat.String

if VERSION < v"0.5.0-dev+2915"
    import Compat.issymmetric
else
    import Base.issymmetric
end

using Base.LinAlg

import Base.show

export PardisoSolver, MKLPardisoSolver
export set_iparm!, set_dparm!, set_matrixtype!, set_solver!, set_phase!, set_msglvl!, set_nprocs!
export get_iparm, get_iparms, get_dparm, get_dparms
export get_matrixtype, get_solver, get_phase, get_msglvl, get_nprocs
export set_maxfct!, set_perm!, set_mnum!
export get_maxfct, get_perm, get_mnum
export checkmatrix, checkvec, printstats, pardisoinit, pardiso
export solve, solve!
export get_matrix

macro R_str(s)
    s
end

include("CEnum.jl")

using .CEnum

type PardisoException <: Exception
    info::String
end
type PardisoPosDefException <: Exception
    info::String
end
Base.showerror(io::IO, e::Union{PardisoException, PardisoPosDefException}) = print(io, e.info);

typealias PardisoNumTypes Union{Float64, Complex128}

abstract AbstractPardisoSolver

function __init__()
    # Global variables used here are defined in the created deps.jl file in the deps folder
    if !(MKL_PARDISO_LIB_FOUND || PARDISO_LIB_FOUND)
        warn(string("""No Pardiso library found when Pkg.build("Pardiso") ran, this package will not currently be usable. """,
                    """Please see the installation instructions and rerun Pkg.build("Pardiso")."""))
    end

   if MKL_PARDISO_LIB_FOUND
        try
            @static if is_windows() begin
                global const libmkl_core = Libdl.dlopen(joinpath(MKLROOT, "..", "redist", "intel64", "mkl", "mkl_rt.dll"), Libdl.RTLD_GLOBAL)
                global const mkl_init = Libdl.dlsym(libmkl_core, "pardisoinit")
                global const mkl_pardiso_f = Libdl.dlsym(libmkl_core, "pardiso")
                global const set_nthreads = Libdl.dlsym(libmkl_core, "mkl_domain_set_num_threads")
                global const get_nthreads = Libdl.dlsym(libmkl_core, "mkl_domain_get_max_threads")
            end
            end

            # Using the libmkl_rt on Ubuntu hung my computer
            @static if is_unix() begin
                global const libgomp = Libdl.dlopen("libgomp", Libdl.RTLD_GLOBAL)
                global const libmkl_core = Libdl.dlopen(string(MKLROOT, "/lib/intel64/libmkl_core"), Libdl.RTLD_GLOBAL)
                global const libmkl_threaded = Libdl.dlopen(string(MKLROOT, "/lib/intel64/libmkl_gnu_thread"), Libdl.RTLD_GLOBAL)
                global const libmkl_gd = Libdl.dlopen(string(MKLROOT, "/lib/intel64/libmkl_gf_lp64"), Libdl.RTLD_GLOBAL)
                global const mkl_init = Libdl.dlsym(libmkl_gd, "pardisoinit")
                global const mkl_pardiso_f = Libdl.dlsym(libmkl_gd, "pardiso")
                global const set_nthreads = Libdl.dlsym(libmkl_gd, "mkl_domain_set_num_threads")
                global const get_nthreads = Libdl.dlsym(libmkl_gd, "mkl_domain_get_max_threads")
            end
            end
            global const MKL_PARDISO_LOADED = true

        catch e
            warn("MKL Pardiso did not manage to load, error thrown was: $e")
            global const MKL_PARDISO_LOADED = false
        end
    else
        global const MKL_PARDISO_LOADED = false
    end

    if PARDISO_LIB_FOUND
        try
            global const libpardiso = Libdl.dlopen(PARDISO_PATH, Libdl.RTLD_GLOBAL)
            global const init = Libdl.dlsym(libpardiso, "pardisoinit")
            global const pardiso_f = Libdl.dlsym(libpardiso, "pardiso")
            global const pardiso_chkmatrix = Libdl.dlsym(libpardiso, "pardiso_chkmatrix")
            global const pardiso_chkmatrix_z = Libdl.dlsym(libpardiso, "pardiso_chkmatrix_z")
            global const pardiso_printstats = Libdl.dlsym(libpardiso, "pardiso_printstats")
            global const pardiso_printstats_z = Libdl.dlsym(libpardiso, "pardiso_printstats_z")
            global const pardiso_chkvec = Libdl.dlsym(libpardiso, "pardiso_chkvec")
            global const pardiso_chkvec_z = Libdl.dlsym(libpardiso, "pardiso_chkvec_z")

            # Windows Pardiso lib comes with BLAS + LAPACK prebaked but not on UNIX so we open them here
            # if not MKL is loaded
            @static if is_unix() begin
                if !isdefined(Pardiso, :libgomp)
                    global const libgomp = Libdl.dlopen("libgomp", Libdl.RTLD_GLOBAL)
                end
                if !MKL_PARDISO_LOADED
                    global const libblas = Libdl.dlopen("libblas", Libdl.RTLD_GLOBAL)
                    global const liblapack = Libdl.dlopen("liblapack", Libdl.RTLD_GLOBAL)
                end
                global const libgfortran = Libdl.dlopen("libgfortran", Libdl.RTLD_GLOBAL)
            end
            end
            global const PARDISO_LOADED = true
        catch e
            warn("Pardiso did not manage to load, error thrown was: $e")
            global const PARDISO_LOADED = false
        end
    else
        global const PARDISO_LOADED = false
    end

end

include("../deps/deps.jl")
include("enums.jl")
include("project_pardiso.jl")
include("mkl_pardiso.jl")

# Getters and setters
set_matrixtype!(ps::AbstractPardisoSolver, v::Int) = set_matrixtype!(ps, MatrixType[v][1])
function set_matrixtype!(ps::AbstractPardisoSolver, v::MatrixType)
    ps.mtype = v
end

get_matrixtype(ps::AbstractPardisoSolver) = ps.mtype
get_iparm(ps::AbstractPardisoSolver, i::Integer) = ps.iparm[i]
get_iparms(ps::AbstractPardisoSolver) = ps.iparm
set_iparm!(ps::AbstractPardisoSolver, i::Integer, v::Integer) = ps.iparm[i] = v

get_mnum(ps::AbstractPardisoSolver) = ps.mnum
set_mnum!(ps::AbstractPardisoSolver, mnum::Integer) = ps.mnum = mnum

get_maxfct(ps::AbstractPardisoSolver) = ps.maxfct
set_maxfct!(ps::AbstractPardisoSolver, maxfct::Integer) = ps.maxfct = maxfct

get_perm(ps::AbstractPardisoSolver) = ps.perm
set_perm!{T <: Integer}(ps::AbstractPardisoSolver, perm::Vector{T}) = ps.perm = convert(Vector{Int32}, perm)

get_phase(ps::AbstractPardisoSolver) = ps.phase

set_phase!(ps::AbstractPardisoSolver, v::Int) = set_phase!(ps, Phase[v][1])
function set_phase!(ps::AbstractPardisoSolver, v::Phase)
    ps.phase = v
end

get_msglvl(ps::AbstractPardisoSolver) = ps.msglvl

set_msglvl!(ps::AbstractPardisoSolver, v::Integer) = set_msglvl!(ps, MessageLevel[v][1])
function set_msglvl!(ps::AbstractPardisoSolver, v::MessageLevel)
    ps.msglvl = v
end

function pardisoinit(ps::AbstractPardisoSolver)
    ccall_pardisoinit(ps)
    return
end

function solve{Ti, Tv <: PardisoNumTypes}(ps::AbstractPardisoSolver, A::SparseMatrixCSC{Tv, Ti},
                                       B::VecOrMat{Tv}, T::Symbol=:N)
  X = copy(B)
  solve!(ps, X, A, B, T)
  return X
end

function solve!{Ti, Tv <: PardisoNumTypes}(ps::AbstractPardisoSolver, X::VecOrMat{Tv},
                                        A::SparseMatrixCSC{Tv, Ti}, B::VecOrMat{Tv},
                                        T::Symbol=:N)

    pardisoinit(ps)

    # We need to set the transpose flag in PARDISO when we DON'T want
    # a transpose in Julia because we are passing a CSC formatted
    # matrix to PARDISO which expects a CSR matrix.
    if T == :N
        if isa(ps, PardisoSolver)
            set_iparm!(ps, 12, 1)
        else
            set_iparm!(ps, 12, 2)
        end
    elseif T == :C || T == :T
        set_iparm!(ps, 12, 0)
    else
        throw(ArgumentError("only :T, :N and :C, are valid transpose symbols"))
    end

    # If we want the matrix to only be transposed and not conjugated
    # we have to conjugate it before sending it to Pardiso due to CSC CSR
    # mismatch.

    # This is the heuristics for choosing what matrix type to use
    ##################################################################
    # - If hermitian try to solve with symmetric positive definite.
    #   - On pos def exception, solve instead with symmetric indefinite.
    # - If complex and symmetric, solve with symmetric complex solver
    # - Else solve as unsymmetric.
    if ishermitian(A)
        eltype(A) == Float64 ? set_matrixtype!(ps, REAL_SYM_POSDEF) : set_matrixtype!(ps, COMPLEX_HERM_POSDEF)
        try
            if typeof(ps) == PardisoSolver
                   pardiso(ps, X, get_matrix(ps, A, T), B)
            else
                # Workaround for #3
                if eltype(A) == Complex128
                    throw(PardisoPosDefException(""))
                end
                pardiso(ps, X, get_matrix(ps, A, T), B)
            end
        catch e
            isa(e, PardisoPosDefException) || rethrow(e)
            eltype(A) == Float64 ? set_matrixtype!(ps, REAL_SYM_INDEF) : set_matrixtype!(ps, COMPLEX_HERM_INDEF )
            pardiso(ps, X, get_matrix(ps, A, T), B)
        end
    elseif issymmetric(A)
        set_matrixtype!(ps, COMPLEX_SYM)
        pardiso(ps, X, get_matrix(ps, A, T), B)
    else
        eltype(A) == Float64 ? set_matrixtype!(ps, REAL_NONSYM) : set_matrixtype!(ps, COMPLEX_NONSYM)
        pardiso(ps, X, get_matrix(ps, A, T), B)
    end
    original_phase = get_phase(ps)

    # Release memory, TODO: We are running the convert on IA and JA here
    # again which is unnecessary.
    set_phase!(ps, RELEASE_ALL)
    pardiso(ps, X, A, B)
    set_phase!(ps, original_phase)
    return X
end

function pardiso{Ti, Tv <: PardisoNumTypes}(ps::AbstractPardisoSolver, X::VecOrMat{Tv}, A::SparseMatrixCSC{Tv, Ti},
                                            B::VecOrMat{Tv})
    if length(X) != 0
        dim_check(X, A, B)
    end

    if Tv <: Complex && isreal(get_matrixtype(ps))
        throw(ErrorException(string("input matrix is complex while PardisoSolver ",
                                    "has a real matrix type set: $(get_matrixtype(ps))")))
    end

    if Tv <: Real && !isreal(get_matrixtype(ps))
        throw(ErrorException(string("input matrix is real while PardisoSolver ",
                                    "has a complex matrix type set: $(get_matrixtype(ps))")))
    end

    N = Int32(size(A, 2))

    AA = A.nzval
    IA = convert(Vector{Int32}, A.colptr)
    JA = convert(Vector{Int32}, A.rowval)

    NRHS = Int32(size(B, 2))

    ccall_pardiso(ps, N, AA, IA, JA, NRHS, B, X)
end

pardiso(ps::AbstractPardisoSolver) = ccall_pardiso(ps, 0, Float64[], Int32[], Int32[], 0, Float64[], Float64[])
function pardiso{Ti, Tv <: PardisoNumTypes}(ps::AbstractPardisoSolver, A::SparseMatrixCSC{Tv, Ti}, B::VecOrMat{Tv})
    pardiso(ps, Tv[], A, B)
end

function dim_check(X, A, B)
    size(X) == size(B) || throw(DimensionMismatch(string("solution has $(size(X)), ",
                                                         "RHS has size as $(size(B)).")))
    size(A, 1) == size(B, 1) || throw(DimensionMismatch(string("matrix has $(size(A,1)) ",
                                                               "rows, RHS has $(size(B,1)) rows.")))
end

end # module


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    nodedegree,
    sssr, sssrmembership, sssrbondmembership,
    fusedrings, fusedringmembership,
    smallestsssr, sssrcount,
    isringatom, isringbond,
    atomsymbol, charge, multiplicity, bondorder,
    valence, lonepair,
    heavyatomconnected, explicithconnected, implicithconnected,
    hydrogenconnected, connectivity,
    ishdonor, hacceptorcount, ishacceptor, hdonorcount,
    isrotatable, rotatablecount,
    atomcounter, heavyatomcount, molecularformula, empiricalformula,
    pielectron, hybridization,
    isaromaticring, isaromatic, isaromaticbond,
    precalculate!

const LONEPAIR_COUNT = Dict(
    :H => 0, :B => -1, :C => 0, :N => 1, :O => 2, :F => 3,
    :Si => 0, :P => 1, :S => 2, :Cl => 3,
    :As => 1, :Se => 2, :Br => 3, :I => 3
)

# Molecular graph topology

"""
    nodedegree(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing the node degree of the molecular graph
of 1 to ``n``th atoms of the given molecule.

This property corresponds to SMARTS `D` query.
"""
@cachefirst nodedegree(mol::GraphMol) = [degree(mol, n) for n in 1:nodecount(mol)]
nodedegree(view::SubgraphView) = nodedegree(view.graph)

"""
    sssr(mol::GraphMol) -> Vector{Vector{Int}}

Return vectors of ring nodes representing small set of smallest rings (SSSR).

See [`Graph.minimumcyclebasis`](@ref).
"""
@cachefirst function sssr(mol::GraphMol)
    cycles = Vector{Int}[]
    for cy in minimumcyclebasis(mol)
        cp = copy(cy)
        (u, v) = getedge(mol, pop!(cp))
        push!(cycles, shortestpathnodes(edgesubgraph(mol, cp), u, v))
    end
    return cycles
end

"""
    sssrmembership(mol::GraphMol) -> Vector{Set{Int}}

Return a vector of size ``n`` representing [`sssr`](@ref) membership of
1 to ``n``th atoms of the given molecule.

SSSR membership is represented as a set of SSSR indices assigned to each rings.
This means atoms that have the same SSSR index belong to the same SSSR.
"""
@cachefirst function sssrmembership(mol::GraphMol)
    nodes = [Set{Int}() for n in 1:nodecount(mol)]
    for (i, cyc) in enumerate(sssr(mol))
        for n in cyc
            push!(nodes[n], i)
        end
    end
    return nodes
end

"""
    sssrbondmembership(mol::GraphMol) -> Vector{Set{Int}}

Return a vector of size ``n`` representing [`sssr`](@ref) membership of
1 to ``n``th bonds of the given molecule.

SSSR membership is represented as a set of SSSR indices assigned to each rings.
This means bonds that have the same SSSR index belong to the same SSSR.
"""
@cachefirst function sssrbondmembership(mol::GraphMol)
    edges = [Set{Int}() for n in 1:edgecount(mol)]
    for (i, cyc) in enumerate(minimumcyclebasis(mol))
        for e in cyc
            push!(edges[e], i)
        end
    end
    return edges
end

"""
    fusedrings(mol::UndirectedGraph) -> Vector{Set{Int}}

Return vectors of fused ring node sets.

A fused ring is defined as a 2-edge connected components in terms of graph theory.
Spirocyclic structures are considered to be part of a fused ring.
"""
@cachefirst function fusedrings(mol::GraphMol)
    cobr = setdiff(edgeset(mol), bridges(mol))
    subg = plaingraph(edgesubgraph(mol, cobr))
    return connectedcomponents(subg)
end

"""
    fusedringmembership(mol::UndirectedGraph) -> Vector{Int}

Return a vector of size ``n`` representing [`fusedrings`](@ref) membership of
1 to ``n``th atoms of the given molecule.

Fused ring membership is represented as a set of fused ring indices assigned to each fused rings.
This means atoms that have the same fused ring index belong to the same fused ring.
"""
@cachefirst function fusedringmembership(mol::GraphMol)
    arr = [Set{Int}() for i in 1:nodecount(mol)]
    for (i, conn) in enumerate(fusedrings(mol))
        for n in conn
            push!(arr[n], i)
        end
    end
    return arr
end

"""
    smallestsssr(mol::UndirectedGraph) -> Vector{Int}

Return a vector of size ``n`` representing the size of the smallest [`sssr`](@ref)
that 1 to ``n``th atoms of the given molecule belong to.

This property corresponds to SMARTS `r` query.
"""
@cachefirst function smallestsssr(mol::GraphMol)
    sssr_ = sssr(mol)
    arr = []
    for ks in sssrmembership(mol)
        val = isempty(ks) ? 0 : minimum([length(sssr_[k]) for k in ks])
        push!(arr, val)
    end
    return arr
end
smallestsssr(view::SubgraphView) = smallestsssr(view.graph)

"""
    sssrcount(mol::UndirectedGraph) -> Vector{Int}

Return a vector of size ``n`` representing the number of [`sssr`](@ref)
that 1 to ``n``th atoms of the given molecule belong to.

This property corresponds to SMARTS `R` query.
"""
@cachefirst sssrcount(mol::UndirectedGraph) = length.(sssrmembership(mol))
sssrcount(view::SubgraphView) = sssrcount(view.graph)

"""
    isringatom(mol::UndirectedGraph) -> Vector{Bool}

Return a vector of size ``n`` representing whether 1 to ``n``th atoms of
the given molecule belong to a ring or not.
"""
@cachefirst isringatom(mol::UndirectedGraph) = .!isempty.(sssrmembership(mol))
isringatom(view::SubgraphView) = isringatom(view.graph)

"""
    isringbond(mol::UndirectedGraph) -> Vector{Bool}

Return a vector of size ``n`` representing whether 1 to ``n``th bonds of
the given molecule belong to a ring or not.
"""
@cachefirst isringbond(mol::UndirectedGraph) = .!isempty.(sssrbondmembership(mol))
isringbond(view::SubgraphView) = isringbond(view.graph)

# Elemental properties

"""
    atomsymbol(mol::GraphMol) -> Vector{Symbol}

Return a vector of size ``n`` representing atom symbols of 1 to ``n``th atoms of
the given molecule.
"""
@cachefirst atomsymbol(mol::UndirectedGraph) = getproperty.(nodeattrs(mol), :symbol)

"""
    charge(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing atom charges of 1 to ``n``th atoms of
the given molecule.
"""
@cachefirst charge(mol::UndirectedGraph) = getproperty.(nodeattrs(mol), :charge)

"""
    multiplicity(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing atom multiplicities of 1 to ``n``th atoms of
the given molecule (1: non-radical, 2: radical, 3: biradical).
"""
@cachefirst multiplicity(mol::UndirectedGraph) = getproperty.(nodeattrs(mol), :multiplicity)

"""
    bondorder(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing bond order of 1 to ``n``th bonds of
the given molecule.
"""
@cachefirst bondorder(mol::UndirectedGraph) = getproperty.(edgeattrs(mol), :order)

# Valence

"""
    lonepair(mol::GraphMol) -> Vector{Union{Int,Nothing}}

Return a vector of size ``n`` representing the number of lone pairs of
1 to ``n``th atoms of the given molecule.

The number of lone pair in inorganic atoms would be `nothing`.
The result can take negative value if the atom has empty shells (e.g. B).
"""
@cachefirst function lonepair(mol::GraphMol)
    vec = Union{Int,Nothing}[]
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    for i in 1:nodecount(mol)
        num = get(LONEPAIR_COUNT, atomsymbol_[i], nothing)
        v = num === nothing ? nothing : num - charge_[i]
        push!(vec, v)
    end
    return vec
end
lonepair(view::SubgraphView) = lonepair(view.graph)

"""
    apparentvalence(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing the number of total bond order
incident to 1 to ``n``th atoms of the given molecule.
"""
@cachefirst function apparentvalence(mol::GraphMol)
    vec = zeros(Int, nodecount(mol))
    bondorder_ = bondorder(mol)
    for i in 1:nodecount(mol)
        for (inc, adj) in neighbors(mol, i)
            vec[i] += bondorder_[inc]
        end
    end
    return vec
end

apparentvalence(view::SubgraphView) = apparentvalence(view.graph)

"""
    valence(mol::GraphMol) -> Vector{Union{Int,Nothing}}

Return a vector of size ``n`` representing the intrinsic valence of
1 to ``n``th atoms of the given molecule.

The number of implicit hydrogens would be calculated based on the valence.
The valence value in inorganic atoms would be `nothing`.
This property corresponds to SMARTS `v` query.
"""
@cachefirst function valence(mol::GraphMol)
    vec = Union{Int,Nothing}[]
    atomsymbol_ = atomsymbol(mol)
    lonepair_ = lonepair(mol)
    for i in 1:nodecount(mol)
        if lonepair_[i] === nothing
            push!(vec, nothing)
        elseif atomsymbol_[i] === :H
            push!(vec, 1)
        else
            push!(vec, 4 - abs(lonepair_[i]))
        end
    end
    return vec
end
valence(view::SubgraphView) = valence(view.graph)

"""
    explicithconnected(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing the number of explicit hydrogens
connected to 1 to ``n``th atoms of the given molecule.

"Explicit" means hydrogens are connected to the heavy atom as atom nodes.
"""
@cachefirst function explicithconnected(mol::GraphMol)
    vec = zeros(Int, nodecount(mol))
    atomsymbol_ = atomsymbol(mol)
    for i in 1:nodecount(mol)
        for (inc, adj) in neighbors(mol, i)
            if atomsymbol_[adj] === :H
                vec[i] += 1
            end
        end
    end
    return vec
end
explicithconnected(view::SubgraphView) = explicithconnected(view.graph)

"""
    implicithconnected(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing the number of implicit hydrogens
connected to 1 to ``n``th atoms of the given molecule.

"Implicit" means hydrogens are not represented as graph nodes,
but it can be infered from the intrinsic valence of typical organic atoms.
"""
@cachefirst function implicithconnected(mol::UndirectedGraph)
    hcnt = (v, av) -> v === nothing ? 0 : max(0, v - av)
    return hcnt.(valence(mol), apparentvalence(mol))
end

"""
    heavyatomconnected(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing the number of non-hydrogen atoms
connected to 1 to ``n``th atoms of the given molecule.
"""
@cachefirst heavyatomconnected(mol::UndirectedGraph) = nodedegree(mol) - explicithconnected(mol)

"""
    hydrogenconnected(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing the number of total hydrogens
(implicit and explicit) connected to 1 to ``n``th atoms of the given molecule.

This property corresponds to SMARTS `H` query.
"""
@cachefirst hydrogenconnected(mol::UndirectedGraph) = explicithconnected(mol) + implicithconnected(mol)

"""
    connectivity(mol::GraphMol) -> Vector{Int}

Return a vector of size ``n`` representing the number of total atoms
(implicit and explicit) connected to 1 to ``n``th atoms of the given molecule.

This property corresponds to SMARTS `X` query.
"""
@cachefirst connectivity(mol::UndirectedGraph) = nodedegree(mol) + implicithconnected(mol)

# Hydrogen bond donor/acceptor

@cachefirst function ishacceptor(mol::UndirectedGraph)
    ac = (sym, lp) -> lp === nothing ? false : sym in (:N, :O, :F) && lp > 0
    return ac.(atomsymbol(mol), lonepair(mol))
end

"""
    hacceptorcount(mol::GraphMol) -> Int

Return the total number of hydrogen bond acceptors (N, O and F).
"""
hacceptorcount(mol::GraphMol) = reduce(+, ishacceptor(mol); init=0)

@cachefirst function ishdonor(mol::UndirectedGraph)
    dc = (sym, h) -> sym in (:N, :O) && h > 0
    return dc.(atomsymbol(mol), hydrogenconnected(mol))
end

"""
    hdonorcount(mol::GraphMol) -> Int

Return the total number of hydrogen bond donors (O and N attached to hydrogens).
"""
hdonorcount(mol::GraphMol) = reduce(+, ishdonor(mol); init=0)

# Rotatable bonds

"""
    isrotatable(mol::GraphMol)

Return a vector of size ``n`` representing whether 1 to ``n``th bonds
of the given molecule are rotatable or not.
"""
@cachefirst function isrotatable(mol::GraphMol)
    nodedegree_ = nodedegree(mol)
    isringbond_ = isringbond(mol)
    bondorder_ = bondorder(mol)
    vec = Bool[]
    for (i, (u, v)) in enumerate(edgesiter(mol))
        rot = (!isringbond_[i] && bondorder_[i] == 1
            && nodedegree_[u] != 1 && nodedegree_[v] != 1)
        push!(vec, rot)
    end
    return vec
end
isrotatable(view::SubgraphView) = isrotatable(view.graph)

"""
    rotatablecount(mol::GraphMol) -> Int

Return the total number of rotatable bonds.
"""
rotatablecount(mol::GraphMol) = reduce(+, isrotatable(mol); init=0)

# Composition

"""
    atomcounter(mol::GraphMol) -> Dict{Symbol,Int}

Count the number of atoms and return symbol => count dict.
"""
@cachefirst function atomcounter(mol::GraphMol)
    counter = Dict{Symbol,Int}()
    for sym in atomsymbol(mol)
        if !haskey(counter, sym)
            counter[sym] = 1
        else
            counter[sym] += 1
        end
    end
    hcnt = reduce(+, implicithconnected(mol); init=0)
    if hcnt > 0
        if !haskey(counter, :H)
            counter[:H] = hcnt
        else
            counter[:H] += hcnt
        end
    end
    return counter
end

"""
    heavyatomcount(mol::GraphMol) -> Int

Return the total number of non-hydrogen atoms.
"""
function heavyatomcount(mol::GraphMol)
    counter = atomcounter(mol)
    delete!(counter, :H)
    return reduce(+, values(counter); init=0)
end

function writeformula(counter::Dict{Symbol,Int})
    strs = []
    if haskey(counter, :C)
        push!(strs, "C")
        c = pop!(counter, :C)
        c > 1 && push!(strs, string(c))
        if haskey(counter, :H)
            push!(strs, "H")
            h = pop!(counter, :H)
            h > 1 && push!(strs, string(h))
        end
    end
    for sym in sort(collect(keys(counter)))
        push!(strs, string(sym))
        cnt = counter[sym]
        cnt > 1 && push!(strs, string(cnt))
    end
    return join(strs)
end

"""
    molecularformula(mol::GraphMol) -> String

Return the molecular formula in Hill system.
"""
function molecularformula(mol::GraphMol)
    counter = atomcounter(mol)
    return writeformula(counter)
end

"""
    empiricalformula(mol::GraphMol) -> String

Return the empirical formula in Hill system.
"""
function empiricalformula(mol::GraphMol)
    counter = atomcounter(mol)
    if length(counter) > 1
        dv = reduce(gcd, values(counter))
        for k in keys(counter)
            counter[k] = div(counter[k], dv)
        end
    end
    return writeformula(counter)
end

# Hybridization

"""
    pielectron(mol::GraphMol) -> Vector{Int}

Returns a vector of size ``n`` representing the number of ``\\pi`` electrons
of 1 to ``n``th atoms of the given molecule.

The counting of ``\\pi`` electrons is based on the following rules.

- Any atom incident to a double bond -> +1
- Any atom incident to two double bond -> +2
- Any atom incident to a triple bond -> +2
- Any other uncharged N, O or S that are neighbor of multiple bonds -> +2

These rules are applied for only typical organic atoms. The values for inorganic atoms will be 0.
"""
@cachefirst function pielectron(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    pie_ = apparentvalence(mol) - nodedegree(mol)
    vec = zeros(Int, nodecount(mol))
    for i in 1:nodecount(mol)
        vec[i] = pie_[i]
        atomsymbol_[i] in (:N, :O, :S) || continue
        for (inc, adj) in neighbors(mol, i)
            if pie_[i] == 0 && pie_[adj] > 0 && charge_[i] == 0
                vec[i] = 2
                break
            end
        end
    end
    return vec
end
pielectron(view::SubgraphView) = pielectron(view.graph)

"""
    hybridization(mol::GraphMol) -> Vector{Int}

Returns a vector of size ``n`` representing the orbital hybridization symbols
(`:sp3`, `:sp2`, `:sp` or `:none`) of 1 to ``n``th atoms of the given molecule.

The hybridization value in inorganic atoms and non-typical organic atoms will be `:none`
(e.g. s, sp3d and sp3d2 orbitals).
"""
@cachefirst function hybridization(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    pielectron_ = pielectron(mol)
    connectivity_ = connectivity(mol)
    lonepair_ = lonepair(mol)
    vec = Symbol[]
    for i in 1:nodecount(mol)
        if lonepair_[i] === nothing || atomsymbol_[i] === :H
            push!(vec, :none)
            continue
        end
        orbitals = connectivity_[i] + lonepair_[i]
        if orbitals == 4
            if (atomsymbol_[i] in [:N, :O] && pielectron_[i] == 2)
                push!(vec, :sp2)  # adjacent to conjugated bonds
            else
                push!(vec, :sp3)
            end
        elseif orbitals == 3
            push!(vec, :sp2)
        elseif orbitals == 2
            push!(vec, :sp)
        else
            push!(vec, :none)
        end
    end
    return vec
end
hybridization(view::SubgraphView) = hybridization(view.graph)

# Aromaticity

@cachefirst function isaromaticring(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    nodedegree_ = nodedegree(mol)
    pie_ = apparentvalence(mol) - nodedegree(mol)
    lonepair_ = lonepair(mol)
    carbonylO = findall(
        (atomsymbol_ .=== :O) .* (nodedegree_ .== 1) .* (pie_ .== 1))
    carbonylC = falses(nodecount(mol))
    for o in carbonylO
        c = pop!(adjacencies(mol, o))
        if atomsymbol_[c] === :C
            carbonylC[c] = true
        end
    end
    arr = falses(circuitrank(mol))
    for (i, ring) in enumerate(sssr(mol))
        cnt = 0
        for r in ring
            if carbonylC[r]
                continue
            elseif pie_[r] == 1
                cnt += 1
            elseif lonepair_[r] === nothing
                cnt = 0
                break
            elseif lonepair_[r] > 0
                cnt += 2
            elseif lonepair_[r] < 0
                continue
            else
                cnt = 0
                break
            end
        end
        if cnt % 4 == 2
            # Huckel rule check
            arr[i] = true
            continue
        end
        if nodeattrtype(mol) === SmilesAtom
            # SMILES aromatic atom
            sub = nodesubgraph(mol, Set(ring))
            if all(nodeattr(mol, n).isaromatic for n in nodeset(sub))
                arr[i] = true
            end
        end
    end
    return arr
end

"""
    isaromatic(mol::GraphMol) -> Vector{Bool}

Returns a vector of size ``n`` representing whether 1 to ``n``th atoms
of the given molecule belong to an aromatic ring or not.

Some kind of aromaticity resulting from long conjugated chains and charge
delocalization may be unrecognizable. Also, non-classical aromaticity
such as Moebius aromaticity is not considered.
"""
@cachefirst function isaromatic(mol::GraphMol)
    aromatic = falses(nodecount(mol))
    for ring in sssr(mol)[isaromaticring(mol)]
        sub = nodesubgraph(mol, Set(ring))
        for n in nodeset(sub)
            aromatic[n] = true
        end
    end
    return aromatic
end
isaromatic(view::SubgraphView) = isaromatic(view.graph)

"""
    isaromaticbond(mol::GraphMol) -> Vector{Bool}

Returns a vector of size ``n`` representing whether 1 to ``n``th bonds
of the given molecule belong to an aromatic ring or not.

See [`isaromatic`](@ref).
"""
@cachefirst function isaromaticbond(mol::GraphMol)
    aromaticbond = falses(edgecount(mol))
    for ring in sssr(mol)[isaromaticring(mol)]
        sub = nodesubgraph(mol, Set(ring))
        for e in edgeset(sub)
            aromaticbond[e] = true
        end
    end
    return aromaticbond
end
isaromaticbond(view::SubgraphView) = isaromaticbond(view.graph)

"""
    precalculate!(mol::GraphMol)

Convenient method to pre-calculate and cache performance bottleneck descriptors.
"""
function precalculate!(mol)
    setcache!(mol, :minimumcyclebasis)
    setcache!(mol, :sssr)
    setcache!(mol, :lonepair)
    setcache!(mol, :apparentvalence)
    setcache!(mol, :valence)
    setcache!(mol, :isaromaticring)
    nodeattrtype(mol) === SmilesAtom && setcache!(mol, :coordgen)
end


text:
# ATTENTION: This file was generated, do not edit by hand.
module GAP_pkg_unipot

using Pkg.Artifacts
using GAP

function __init__()
    # ensure GAP will load the correct instance of this package
    path = artifact"unipot"
    @debug "GAP package 'unipot' path = " * path
    GAP.Globals.SetPackagePath(GapObj("unipot"), GapObj(path))
end

function load()
    res = GAP.Globals.LoadPackage(GapObj("unipot"))
    if res != true
        error("failed to load GAP package unipot")
    end
end

end # module


text:
"""
This can be used as a CellField as long as one evaluates it
on the stored CellPoint.
"""
struct CellState{T,P<:CellPoint} <: CellField
  points::P
  values::AbstractArray

  function CellState{T}(::UndefInitializer,points::CellPoint) where T
    values = _init_values(T,get_data(points))
    P = typeof(points)
    new{T,P}(points,values)
  end

  function CellState(v::Number,points::CellPoint)
    values = _init_values(v,get_data(points))
    T = typeof(v)
    P = typeof(points)
    new{T,P}(points,values)
  end
end

function _init_values(::Type{T},x::AbstractArray{<:Point}) where T
  N = ndims(x)
  Array{T,N}(undef,size(x))
end

function _init_values(::Type{T},x::AbstractArray{<:AbstractVector{<:Point}}) where T
  [Vector{T}(undef,length(xi)) for xi in x]
end

function _init_values(v::Number,x::AbstractArray{<:Point})
  fill(v,size(x))
end

function _init_values(v::Number,x::AbstractArray{<:AbstractVector{<:Point}})
  [fill(v,length(xi)) for xi in x]
end

function get_data(f::CellState)
  @unreachable """\n
  get_data cannot be called on a CellState

  If you see this error messase it is likelly that you are trying to perform
  an operation that does not make sense for a CellState. In most cases,
  to do the wanted operation, you would need first to project the CellState
  to a FESpace (e.g. via a L2 projection).
  """
end

get_triangulation(f::CellState) = get_triangulation(f.points)
DomainStyle(::Type{CellState{T,P}}) where {T,P} = DomainStyle(P)

function evaluate!(cache,f::CellState,x::CellPoint)
  if f.points === x
    f.values
  else
    @unreachable """\n
    It is not possible to evaluate the given CellState on the given CellPoint.

    a CellState can only be evaluated at the CellPoint it was created from.
    If you want to evaluate at another location, you would need first to project the CellState
    to a FESpace (e.g. via a L2 projection).
    """
  end
end

function CellState{T}(::UndefInitializer,a) where T
  points = get_cell_points(a)
  CellState{T}(undef,points)
end

function CellState(v::Number,a)
  points = get_cell_points(a)
  CellState(v,points)
end

function update_state!(updater::Function,f::CellField...)
  ids = findall(map(i->isa(i,CellState),f))
  @assert length(ids) > 0 """\n
  At least one CellState object has to be given to the update_state! function
  """
  a = f[ids]
  x = first(a).points
  @assert all(map(i->i.points===x,a)) """\n
  All the CellState objects given to the update_state! function need to be
  defined on the same CellPoint.
  """
  fx = map(i->evaluate(i,x),f)
  if num_cells(x) > 0
    fxi = map(first,fx)
    fxiq = map(first,fxi)
    need_to_update, states = first_and_tail(updater(fxiq...))
    @assert isa(need_to_update,Bool) && isa(states,Tuple{Vararg{Number}}) """\n
    Wrong return value of the user-defined updater Function. The signature is

        need_to_update, state = first_and_tail(updater(args...))

    where need_to_update is a Bool telling if we need to update_state and
    states is a Tuple of Number objects with the new states.
    See `first_and_tail` for further details.
    """
    msg = """\n
    The number of new states given by the updater Function does not match
    the number of CellState objects given as arguments in the update_state! Funciton.
    """
    @check length(states) <= length(a) msg
  end
  caches = map(array_cache,fx)
  x_data = get_data(x)
  cache_x = array_cache(x_data)
  _update_state_variables!(updater,caches,fx,cache_x,x_data)
  nothing
end

@noinline function  _update_state_variables!(updater,caches,fx,cache_x,x)
  ncells = length(x)
  for cell in 1:ncells
    fxi = map((c,f)->getindex!(c,f,cell),caches,fx)
    xi = getindex!(cache_x,x,cell)
    for q in 1:length(xi)
      fxiq = map(f->f[q],fxi)
      need_to_update, states = first_and_tail(updater(fxiq...))
      if need_to_update
        _update_states!(fxi,q,states,Val{length(states)}())
      end
    end
  end
end

function _update_states!(b,q,states,::Val{i}) where i
  _update_state!(b,q,states,Val{i}())
  _update_states!(b,q,states,Val{i-1}())
  nothing
end

function _update_states!(b,q,states,::Val{0})
  nothing
end

function _update_state!(b,q,states,::Val{i}) where i
  b[i][q] = states[i]
  nothing
end


text:
import QMTK.QuCircuit
using Compat.Test


text:
MusicXMLTypes = [:ScorePartwise, :Part, :Measure, :Note, :Unpitched, :Rest, :Pitch, :Attributes, :Time, :Transpose, :Clef, :Key, :PartList, :ScorePart, :MidiInstrument, :MidiDevice, :ScoreInstrument]

import TreeViews: hastreeview, numberofnodes, treelabel, treenode
import Base: inferencebarrier, show_circular
for T in MusicXMLTypes
    @eval begin
        function Base.show(io::IO, x::$T)
            # code from Base
            t = typeof(x)
            show(io, inferencebarrier(t))
            print(io, '(')
            nf = nfields(x)
            nb = sizeof(x)
            if nf != 0 || nb == 0
                if !show_circular(io, x)
                    recur_io = IOContext(io, Pair{Symbol,Any}(:SHOWN_SET, x),
                                         Pair{Symbol,Any}(:typeinfo, Any))
                    for i in 1:nf-1 # skip last
                        f = fieldname(t, i)
                        if !isdefined(x, f)
                            print(io, undef_ref_str)
                        else
                            show(recur_io, getfield(x, i))
                        end
                        if i < nf
                            print(io, ", ")
                        end
                    end
                    print(io, "aml")
                end
            else
                print(io, "0x")
                r = Ref(x)
                GC.@preserve r begin
                    p = unsafe_convert(Ptr{Cvoid}, r)
                    for i in (nb - 1):-1:0
                        print(io, string(unsafe_load(convert(Ptr{UInt8}, p + i)), base = 16, pad = 2))
                    end
                end
            end
            print(io,')')
        end
        # Juno
        Base.show(io::IO, ::MIME"application/prs.juno.inline", x::$T) = x
        # hastreeview(::$T) = true
        numberofnodes(::$T) = fieldcount($T)-1
    end
end


text:
# Internal searching methods

import AxisArrays: searchsortednearest
@test searchsortednearest(1:5, 2.5) === 3
@test searchsortednearest(1:5, prevfloat(2.5)) === 2
@test searchsortednearest([1,1,2,2,3,3], 1.5) === 3
@test searchsortednearest([1,1,2,2,3,3], 2.0) === 3
@test searchsortednearest([1,1,2,2,3,3], 2.4) === 4
@test searchsortednearest([1,1,2,2,3,3], 3.0) === 5
@test searchsortednearest([1,1,2,2,3,3], Inf) === 6
@test searchsortednearest([1,1,2,2,3,3], -Inf) === 1

# unsafe searchsorted for ranges
import AxisArrays: unsafe_searchsorted
@test unsafe_searchsorted(1:10, -1 .. 1) === -1:1
@test unsafe_searchsorted(1:10, 12 .. 15) === 12:15
@test unsafe_searchsorted(0:2:10, -3 .. -1) === 0:0
@test unsafe_searchsorted(0:2:10, -5 .. 3) === -1:2

@test unsafe_searchsorted(1:2, 4.5 .. 4.5) === 5:4
@test unsafe_searchsorted(1:2, 3.5 .. 3.5) === 4:3
@test unsafe_searchsorted(1:2, 2.5 .. 2.5) === 3:2 === searchsorted(1:2, 2.5 .. 2.5)
@test unsafe_searchsorted(1:2, 1.5 .. 1.5) === 2:1 === searchsorted(1:2, 1.5 .. 1.5)
@test unsafe_searchsorted(1:2, 0.5 .. 0.5) === 1:0 === searchsorted(1:2, 0.5 .. 0.5)
@test unsafe_searchsorted(1:2, -0.5 .. -0.5) === 0:-1
@test unsafe_searchsorted(1:2, -1.5 .. -1.5) === -1:-2

@test unsafe_searchsorted(2:2:4, 0x6 .. 0x6) === 3:3
@test unsafe_searchsorted(2:2:4, 0x5 .. 0x5) === searchsorted(2:2:4, 0x5 .. 0x5) === 3:2
@test unsafe_searchsorted(2:2:4, 0x4 .. 0x4) === searchsorted(2:2:4, 0x4 .. 0x4) === 2:2
@test unsafe_searchsorted(2:2:4, 0x3 .. 0x3) === searchsorted(2:2:4, 0x3 .. 0x3) === 2:1
@test unsafe_searchsorted(2:2:4, 0x2 .. 0x2) === searchsorted(2:2:4, 0x2 .. 0x2) === 1:1
@test unsafe_searchsorted(2:2:4, 0x1 .. 0x1) === searchsorted(2:2:4, 0x1 .. 0x1) === 1:0
@test unsafe_searchsorted(2:2:4, 0x0 .. 0x0) === 0:0