text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# This file contains utility functions

function df_to_dict_with_id(df::DataFrame)
    Dict(col[1] => Dict(zip(df[1], col[2])) for col in eachcol(df, true))
end

function df_to_dict(df::DataFrame)
    Dict(col[1] => col[2] for col in eachcol(df, true))
end

function dictzip(df::DataFrame, which::Tuple{Symbol, Symbol})
    return Dict(zip(df[which[1]], df[which[2]]))
end

function insert_default!(dict::Dict, default, keys)
	for k in keys
		haskey(dict, k) || (dict[k] = default)
	end
end

function jumpvar_to_df(jv::JuMP.JuMPArray;
		dim_names::Array{Symbol, 1} = Symbol[],
		dual::Bool=false)
	dim_arr = jv.indexsets
	dims = length(dim_arr)
	if dual
		arr = getdual(jv).innerArray
	else
		arr = getvalue(jv).innerArray
	end

	length(dim_names) == 0 && (dim_names = [Symbol("x$i") for i in 1:ndims(arr)])

	rows = []
	for ind in CartesianIndices(size(arr))
		row_ind = [dim_arr[dim][ind.I[dim]] for dim in 1:dims]
		push!(rows, (row_ind..., arr[ind]))
	end
	rows = vcat(rows...)

	k = [dim_names[i] => [row[i] for row in rows] for i in 1:length(dim_names)]
	kv = vcat(k..., :Value => [row[length(row)] for row in rows])

	return DataFrame(kv...)
end

get_variable_names(m::JuMP.Model) = [v.name for (k,v) in m.varData]

week_slices(sec_w::Int) = [1:sec_w-1, [t+167 > 8760 ? (t:8760) : (t:t+167) for t in sec_w:168:8760]...]


text:
using OctreeBH
import OctreeBH: AbstractData, assign_additional_node_data!

using StaticArrays

mutable struct Data2{N,T<:Real} <: AbstractData{N,T}
	pos::SVector{N,T}
	idx::Int64
	hsml::T
	mass::T
	#define additional members below
	massH2::T
end

Data2{N,T}() where {N,T} = Data2{N,T}(zero(SVector{N,T}),0,0,0,0)

function assign_additional_node_data!(n::Data2, old::Data2, new::Data2)
	n.massH2 = old.massH2 + new.massH2
end

const BOXSIZE = 1.0

const N = 3 #spatial dimension
const T = Float64

const Npart = 100000 #number of particles

#desired number of neighbor particles (ngbs) within a search radius (hsml)
const Nngb0 = 32

#searching radius that should contain roughly Nngb0 neighboring particles
const hsml0 = BOXSIZE * (Nngb0/(4*pi/3*Npart))^(1/N)

#randomly distributing particles
X = [@SVector rand(N) for _ in 1:Npart]

topnode_length = @SVector(ones(N)) * BOXSIZE  #actual length of tree

center = topnode_length .* 0.5

#individual smoothing length (only relevant for scatter ngbs)
δ = 0.0 #add a small perturbation
hsml = ones(Npart) .* hsml0 .* ( 1 .+ δ .* (rand(Npart) .- 0.5) )

mass = ones(Npart) #particle mass

part = [Data2{N,T}(SVector(X[i]), i, hsml[i], mass[i], 0.6*mass[i]) for i in eachindex(X)]

#build the tree
tree = buildtree(part, center, topnode_length);

#search center
x = @SVector rand(N)

boxsizes = @SVector(ones(N)) * BOXSIZE  #for periodic B.C.

#get indices of the gather neighbors
idx_ngbs_g = get_gather_ngb_tree(x, hsml0, tree, boxsizes)

#get indices of the scatter neighbors
idx_ngbs_s = get_scatter_ngb_tree(x, tree, boxsizes)

#for constant hsml, the two ngbs are identical
@show length(idx_ngbs_g), length(idx_ngbs_s)

ext: jl


text:
using Test;

function run(fname)
    f = open(fname)
    lines = readlines(f)
    close(f)

    mask = ""
    mem_values = Dict()
    for line in lines
        parts = split(line, " = ")
        cmd, val = parts[1], parts[2]
        println("cmd = ", cmd)
        println("val = ", val)
        if cmd == "mask"
            # update the mask
            println("MASK VAL = ", val)
            mask = val
        else
            # update something in mem, applying mask
            mem_idx = parse(Int, cmd[5:length(cmd) - 1])
            num = parse(Int, parts[2])
            println("Mem_idx = ", mem_idx, " ... Num = ", num)

            # convert number to bit representation
            bits = bitstring(num)[29:64]  # want the last 36-bits
            # apply mask
            println("mask = ", mask)
            println("bits = ", bits)
            out = ""
            for (idx, m) in enumerate(mask)
                if m == 'X'
                    out *= bits[idx]
                else
                    out *= m
                end
            end
            println("new  = ", out)
            # store final number in memory
            new_val = parse(Int, out, base=2)
            println("newval = ", new_val)
            mem_values[mem_idx] = new_val
        end
    end
    total = 0
    for k in keys(mem_values)
        v = mem_values[k]
        total += v
    end

    println("[Part 1] Result of running ", fname, " was = ", total)
    return total
end

function run2(fname)
    # Bitmask modifies how we write to memory
    f = open(fname)
    lines = readlines(f)
    close(f)

    mask = ""
    mem_values = Dict()
    for line in lines
        parts = split(line, " = ")
        cmd, val = parts[1], parts[2]
        println("cmd = ", cmd)
        println("val = ", val)
        if cmd == "mask"
            # update the mask
            println("MASK VAL = ", val)
            mask = val
        else
            # update something in mem, applying mask
            mem_idx = parse(Int, cmd[5:length(cmd) - 1])
            num = parse(Int, parts[2])
            println("Mem_idx = ", mem_idx, " ... Num = ", num)

            # convert number to bit representation
            bits = bitstring(mem_idx)[29:64]  # want the last 36-bits
            # apply mask
            println("mask = ", mask)
            println("bits = ", bits)
            out = [""]

            # Generate memory address(es)
            for (idx, m) in enumerate(mask)
                if m == 'X'
                    # add both possible characters to each item
                    out1 = map((x) -> x * '0', out)
                    out2 = map((x) -> x * '1', out)
                    out = [out1..., out2...]
                elseif m == '1'
                    # overwrite with '1'
                    out = map((x) -> x * '1', out)
                elseif m == '0'
                    # don't override
                    out = map((x) -> x * bits[idx], out)
                end
            end
            # println("new  = ", out)

            # store number in memory
            for o in out
                addr = parse(Int, o, base=2)
                # println("Writing [", addr, "] => ", num)
                mem_values[addr] = num
            end
        end
    end

    total = 0
    # println(mem_values)
    for k in keys(mem_values)
        v = mem_values[k]
        total += v
    end

    println("[Part 2] Result of running ", fname, " was = ", total)
    return total
end

# @test run("14ex.txt") == 165
# run("14.txt")

@test run2("14ex.part2.txt") == 208
run2("14.txt")


text:
# This example was adapted from: Coding Challenge #64.1 "Foward Kinematics"
# by Daniel Shiffman
#
# https://thecodingtrain.com/CodingChallenges/064.1-forward-kinematics.html

module FowardKinematics

using Lynx
using Luxor
using Utils

Lynx.init("FowardKinematics", 800, 600)

include("segment.jl")
include("main.jl")

run!(update, setup, await=true)

end # FowardKinematics

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
#=
A diffusion equation with different species coupled through boundary conditions only.
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("coupledbc1d");

useLog("coupledbc1dlog")

# Configuration setup
domain(1)
solverType(CG)
functionSpace(order=2)
timeStepper(EULER_IMPLICIT)

# Mesh
n = 40 # number of elements
mesh(LINEMESH, elsperdim=n)

# Variables and BCs
a = variable("a")
b = variable("b")
testSymbol("v")
boundary(a, 1, NEUMANN, 0)
boundary(b, 1, DIRICHLET, "2*a/(1+a*b)")

# Time interval and initial condition
T = 1;
timeInterval(T)
initial(a, "sin(pi*x)^2")
initial(b, 0)

weakForm([a,b], ["Dt(a*v) + 0.01 * dot(grad(a),grad(v))",
                    "Dt(b*v) + 0.01 * dot(grad(b),grad(v))"])

exportCode("coupledbc1dcode") # uncomment to export generated code to a file
#importCode("coupledbc1dcode") # uncomment to import code from a file

solve([a,b])

finalize_finch()

##### Uncomment below to plot
# x = Finch.grid_data.allnodes[1,:];
# using Plots
# pyplot();
# display(plot([x x], [a.values[:] b.values[:]], marker=:circle))

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# SpinImage.jl: interface to Image objects
export SpinImage, bpp, id, offset, padding, save, timestamp

"""
 Spinnaker SDK Image object
"""
mutable struct SpinImage
  handle::spinImage

  function SpinImage(handle)
    image = new(handle)
    finalizer(_release!, image)
    return image
  end
end

unsafe_convert(::Type{spinImage}, image::SpinImage) = image.handle
unsafe_convert(::Type{Ptr{spinImage}}, image::SpinImage) = pointer_from_objref(image)

function SpinImage()
    @assert spinsys.handle != C_NULL
    himage_ref = Ref(spinImage(C_NULL))
    spinImageCreateEmpty(himage_ref)
    return SpinImage(himage_ref[])
end

# Release handle to image
function _release!(image::SpinImage)
    spinImageDestroy(image)
    image.handle = C_NULL
    return nothing
end

"""
  show(::IO, ::SpinImage)

  Write details of camera to supplied IO.
"""
function show(io::IO, image::SpinImage)
    pixfmt = Ref(spinPixelFormatEnums(0))
    spinImageGetPixelFormat(image, pixfmt)
    write(io, "Spinnaker Image, $(size(image)), $(bpp(image))bpp, $(pixfmt[])")
end

# Return size of image data buffer. Note that this may include metadata.
function _buffersize(image::SpinImage)
    sz = Ref(Csize_t(0))
    spinImageGetBufferSize(image, sz)
    return Int(sz[])
end

# Return a pointer to the data buffer of the input image.
function _bufferptr(image::SpinImage)
    ptr = Ref(Ptr{Cvoid}(0))
    spinImageGetData(image, ptr)
    return Ptr{UInt8}(ptr[])
end

"""
    convert(::SpinImage, ::spinPixelFormatEnums) -> ::SpinImage

    Convert input image to specified pixel format and return new image.
"""
function convert(image::SpinImage, fmt::spinPixelFormatEnums)
    out = SpinImage()
    spinImageConvert(out, fmt, image);
    return out
end

"""
    bpp(::SpinImage) -> Int

    Return bits per pixel of input image.
"""
function bpp(image::SpinImage)
    bpp = Ref(Csize_t(0))
    spinImageGetBitsPerPixel(image, bpp)
    return Int(bpp[])
end

"""
    offset(::SpinImage) -> Tuple{Int}

    Return image offset as tuple (x_offset, y_offset).
"""
function offset(image::SpinImage)
    xoff = Ref(Csize_t(0))
    yoff = Ref(Csize_t(0))
    spinImageGetOffsetX(image, xoff)
    spinImageGetOffsetY(image, yoff)
    return (Int(xoff[]), Int(yoff[]))
end

"""
    padding(::SpinImage) -> Tuple{Int}

    Return image padding as tuple (x_pad, y_pad).
"""
function padding(image::SpinImage)
    xpad = Ref(Csize_t(0))
    ypad = Ref(Csize_t(0))
    spinImageGetPaddingX(image, xpad)
    spinImageGetPaddingY(image, ypad)
    return (Int(xpad[]), Int(ypad[]))
end

"""
    size(::SpinImage) -> Tuple{Int}

    Return image size as tuple (width, height).
"""
function size(image::SpinImage)
    width = Ref(Csize_t(0))
    height = Ref(Csize_t(0))
    spinImageGetWidth(image, width)
    spinImageGetHeight(image, height)
    return (Int(width[]), Int(height[]))
end

"""
    timestamp(::SpinImage) -> Int64

    Return image timestamp from chunk data in nanoseconds since the last timeclock reset (i.e. at camera boot).
"""
function timestamp(image::SpinImage)
    timestamp_ref = Ref(Int64(0))
    spinImageChunkDataGetIntValue(image, "ChunkTimestamp", timestamp_ref)
    return timestamp_ref[]
end

"""
    id(::SpinImage) -> Int64

    Return image id from image chunk data.
"""
function id(image::SpinImage)
    id_ref = Ref(Int64(0))
    spinImageChunkDataGetIntValue(image, "ChunkFrameID", id_ref);
    return id_ref[]
end

"""
    exposure(::SpinImage) -> Float64

    Return exposure time from image chunk data.
"""
function exposure(image::SpinImage)
    hExposure = Ref(Float64(0))
    spinImageChunkDataGetFloatValue(image, "ChunkExposureTime", hExposure);
    return hExposure[]
end

"""
    save(fn::AbstractString, ::SpinImage, ::spinImageFileFormat)

    Save the input image to file `fn` in specified format.
"""
function save(image::SpinImage, fn::AbstractString, fmt::spinImageFileFormat)
    spinImageSave(image, fn, fmt)
end


text:
with_repo_access() do test_repo, path
    context("test read note for object") do
        id = Oid("36060c58702ed4c2a40832c51758d5344201d89a")
        obj = test_repo[id]
        n = notes(obj)
        @test message(n) == "note text\n"
        @test Oid(n) == Oid("94eca2de348d5f672faf56b0decafa5937e3235e")
    end

    context("test read note for object from ref") do
        id = Oid("36060c58702ed4c2a40832c51758d5344201d89a")
        obj = test_repo[id]
        n = notes(obj, "refs/notes/commits")
        @test message(n) == "note text\n"
        @test Oid(n) == Oid("94eca2de348d5f672faf56b0decafa5937e3235e")
    end

    context("test object without note") do
        id = Oid("8496071c1b46c854b31185ea97743be6a8774479")
        obj = test_repo[id]
        @test notes(obj) == nothing
    end

    context("test nothing ref lookup") do
        id = Oid("36060c58702ed4c2a40832c51758d5344201d89a")
        obj = test_repo[id]
        @test notes(obj, "refs/notes/missing") == nothing
    end

    context("test iterate over notes") do
        n = 0
        for (note, oid) in foreach(GitNote, test_repo, "refs/notes/commits")
            n += 1
            @test message(note) == "note text\n"
            @test oid == Oid("36060c58702ed4c2a40832c51758d5344201d89a")
        end
        @test n == 1
    end

    context("test each note iterable") do
        @test isa(foreach(GitNote, test_repo, "refs/notes/commits"), LibGit2.NoteIterator)
    end

    context("test default ref") do
        @test note_default_ref(test_repo) == "refs/notes/commits"
    end
end

with_tmp_repo_access("test create note") do test_repo, path
    sig = Signature("Julia", "julia@julia.com")
    id = Oid("8496071c1b46c854b31185ea97743be6a8774479")
    msg = "This is the note message\n\nThis note is created from Rugged"
    obj = test_repo[id]

    note_id = create_note!(obj, msg, committer=sig, author=sig, ref="refs/notes/test")

    @test note_id == Oid("38c3a690c474d8dcdb13088205a464a60312eec4")
    # note is actually a blob
    blob = test_repo[note_id]
    @test Oid(blob) == note_id
    @test bytestring(blob) == msg
    @test isa(blob, GitBlob)

    n = notes(obj, "refs/notes/test")
    @test Oid(n) == note_id
    @test message(n) == msg
end

with_tmp_repo_access("test create note without signature") do test_repo, path
    testname = "Julia"
    testemail = "julia@julia.com"
    GitConfig(test_repo)["user.name"] = testname
    GitConfig(test_repo)["user.email"] = testemail

    id = Oid("8496071c1b46c854b31185ea97743be6a8774479")
    msg ="This is the note message\n\nThis note is created from Rugged"
    obj = test_repo[id]

    note_id = create_note!(obj, msg, ref="refs/notes/test")

    @test note_id == Oid("38c3a690c474d8dcdb13088205a464a60312eec4")
    note_ref = lookup_ref(test_repo, "refs/notes/test")
    note_commit = test_repo[target(note_ref)]

    @test testname == committer(note_commit).name
    @test testemail == committer(note_commit).email
    @test testname == author(note_commit).name
    @test testemail == author(note_commit).email
end

with_tmp_repo_access("test create note on object with notes raises exception") do test_repo, path
    sig = Signature("Julia", "julia@julia.com")
    id  = Oid("8496071c1b46c854b31185ea97743be6a8774479")
    msg = "This is the note message\n\nThis note is created from Rugged"
    obj = test_repo[id]

    create_note!(obj, msg, committer=sig, author=sig, ref="refs/notes/test")
    @test_throws LibGitError{:Repo,:Exists} create_note!(obj, msg, committer=sig, author=sig, ref="refs/notes/test")
end

with_tmp_repo_access("test overwrite object note") do test_repo, path
    sig = Signature("Julia", "julia@julia.com")
    id = Oid("8496071c1b46c854b31185ea97743be6a8774479")
    msg ="This is the note message\n\nThis note is created from Rugged"
    obj = test_repo[id]

    create_note!(obj, msg, committer=sig, author=sig, ref="refs/notes/test")
    create_note!(obj, "new msg", committer=sig, author=sig, ref="refs/notes/test", force=true)
    note = notes(obj, "refs/notes/test")
    @test message(note) == "new msg"
end

with_tmp_repo_access("test remove note") do test_repo, path
    id = Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    sig = Signature("Julia", "julia@julia.com")
    msg ="This is the note message\n\nThis note is created from Rugged"
    obj = test_repo[id]

    create_note!(obj, msg, author=sig, committer=sig, ref="refs/notes/test")

    @test remove_note!(obj, committer=sig, author=sig, ref="refs/notes/test") == true
    @test notes(obj, "refs/notes/test") == nothing
end

with_tmp_repo_access("test remote without signature") do test_repo, path
    testname = "Julia"
    testemail = "rugged@example.com"
    GitConfig(test_repo)["user.name"] = testname
    GitConfig(test_repo)["user.email"] = testemail
    id = Oid("36060c58702ed4c2a40832c51758d5344201d89a")

    msg ="This is the note message\n\nThis note is created from Rugged"
    obj = test_repo[id]

    create_note!(obj, msg, ref="refs/notes/test")
    create_note!(obj, msg)

    @test remove_note!(obj, ref="refs/notes/test")
    @test remove_note!(obj)
    @test notes(obj, "refs/notes/test") == nothing
    @test notes(obj) == nothing
end

with_tmp_repo_access("test remove missing note") do test_repo, path
    sig = Signature("Julia", "julia@julia.com")
    id = Oid("36060c58702ed4c2a40832c51758d5344201d89a")
    obj = test_repo[id]
    @test remove_note!(obj, committer=sig, author=sig, ref="refs/notes/test") == false
end


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

## Definitions of special intervals and associated functions

## Empty interval:
doc"""`emptyinterval`s are represented as the interval [∞, -∞]; note
that this interval is an exception to the fact that the lower bound is
larger than the upper one."""
emptyinterval{T<:Real}(::Type{T}) = Interval{T}(Inf, -Inf)
emptyinterval{T<:Real}(x::Interval{T}) = emptyinterval(T)
emptyinterval() = emptyinterval(precision(Interval)[1])
const ∅ = emptyinterval(Float64)

isempty(x::Interval) = x.lo == Inf && x.hi == -Inf

const ∞ = Inf

## Entire interval:
doc"""`entireinterval`s represent the whole Real line: [-∞, ∞]."""
entireinterval{T<:Real}(::Type{T}) = Interval{T}(-Inf, Inf)
entireinterval{T<:Real}(x::Interval{T}) = entireinterval(T)
entireinterval() = entireinterval(precision(Interval)[1])

isentire(x::Interval) = x.lo == -Inf && x.hi == Inf
isunbounded(x::Interval) = x.lo == -Inf || x.hi == Inf

# NaI: not-an-interval
doc"""`NaI` not-an-interval: [NaN, NaN]."""
nai{T<:Real}(::Type{T}) = Interval{T}(convert(T, NaN), convert(T, NaN))
nai{T<:Real}(x::Interval{T}) = nai(T)
nai() = nai(precision(Interval)[1])

isnai(x::Interval) = isnan(x.lo) || isnan(x.hi)

isfinite(x::Interval) = isfinite(x.lo) && isfinite(x.hi)
isnan(x::Interval) = isnai(x)

doc"""`isthin(x)` corresponds to `isSingleton`, i.e. it checks if `x` is the set consisting of a single exactly representable float. Thus any float which is not exactly representable does *not* yield a thin interval."""
function isthin(x::Interval)
    # (m = mid(x); m == x.lo || m == x.hi)
    x.lo == x.hi
end

doc"`iscommon(x)` checks if `x` is a **common interval**, i.e. a non-empty, bounded, real interval."
function iscommon(a::Interval)
    (isentire(a) || isempty(a) || isnai(a) || isunbounded(a)) && return false
    true
end

doc"`widen(x)` widens the lowest and highest bounds of `x` to the previous and next representable floating-point numbers, respectively."
widen{T<:AbstractFloat}(x::Interval{T}) = Interval(prevfloat(x.lo), nextfloat(x.hi))


text:
function get_ax2o_a(obj::DSSTTesseralContext)
    return jcall(obj, "getAx2oA", jdouble, ())
end

function get_bo_a(obj::DSSTTesseralContext)
    return jcall(obj, "getBoA", jdouble, ())
end

function get_bo_a_bpo(obj::DSSTTesseralContext)
    return jcall(obj, "getBoABpo", jdouble, ())
end

function get_chi(obj::DSSTTesseralContext)
    return jcall(obj, "getChi", jdouble, ())
end

function get_chi2(obj::DSSTTesseralContext)
    return jcall(obj, "getChi2", jdouble, ())
end

function get_co2ab(obj::DSSTTesseralContext)
    return jcall(obj, "getCo2AB", jdouble, ())
end

function get_e2(obj::DSSTTesseralContext)
    return jcall(obj, "getE2", jdouble, ())
end

function get_max_ecc_pow(obj::DSSTTesseralContext)
    return jcall(obj, "getMaxEccPow", jint, ())
end

function get_mean_motion(obj::DSSTTesseralContext)
    return jcall(obj, "getMeanMotion", jdouble, ())
end

function get_moa(obj::DSSTTesseralContext)
    return jcall(obj, "getMoa", jdouble, ())
end

function get_oo_ab(obj::DSSTTesseralContext)
    return jcall(obj, "getOoAB", jdouble, ())
end

function get_orbit_period(obj::DSSTTesseralContext)
    return jcall(obj, "getOrbitPeriod", jdouble, ())
end

function get_ratio(obj::DSSTTesseralContext)
    return jcall(obj, "getRatio", jdouble, ())
end

function get_res_orders(obj::DSSTTesseralContext)
    return jcall(obj, "getResOrders", List, ())
end

function get_roa(obj::DSSTTesseralContext)
    return jcall(obj, "getRoa", jdouble, ())
end

function get_theta(obj::DSSTTesseralContext)
    return jcall(obj, "getTheta", jdouble, ())
end


text:
module TestOfNFLScrapR

using  DataFrames: DataFrame
using  Test

using  NFLTables: SeasonPart
using  NFLTables.NFLScrapR: getgamedata, getplaydata

@testset "nflscrapR-data df tests" begin
    let season = 2019
        for part in instances(SeasonPart)
            df = getplaydata(season, part)
            @test isa(df, DataFrame)
            df = getgamedata(season, part)
            @test isa(df, DataFrame)
        end
    end
end

end  # module TestOfNFLScrapR


text:
#! format: off
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/jump-dev/ECOS.jl
#############################################################################
# test/direct.jl
# Test the direct interface for the ECOS.jl solver wrapper
#############################################################################

import ECOS
import SparseArrays
using Test

"""
    test_options()

The values below are copied from data.h in the ECOS source code.
"""
function test_options()
    println(ECOS.ECOS_ver())
    n = 223
    m = 220
    p = 114
    l = 201
    ncones = 6
    c = Float64[0, 0, 0, 0.5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    h = Float64[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    q = [3, 3, 3, 3, 3, 4]
    Gjc = [0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220]
    Gir = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219]
    Gpr = Float64[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]
    Ajc = [0, 111, 222, 322, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 539, 540, 541, 564, 565, 566, 589, 590, 591, 618, 619, 620, 632, 633, 634, 635, 636, 637, 638]
    Air = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 111, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 112, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 110, 113, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 3, 9, 15, 27, 53, 54, 57, 66, 67, 71, 74, 87, 92, 96, 98, 100, 101, 2, 8, 12, 19, 29, 30, 34, 37, 39, 42, 43, 44, 47, 50, 51, 56, 59, 70, 73, 79, 82, 85, 94, 102, 103, 0, 4, 7, 13, 14, 16, 17, 21, 23, 28, 31, 33, 36, 38, 60, 64, 78, 81, 88, 90, 91, 95, 99, 104, 105, 1, 6, 10, 20, 22, 26, 32, 35, 41, 46, 48, 49, 52, 55, 58, 65, 68, 69, 75, 76, 77, 80, 83, 84, 86, 93, 97, 106, 107, 5, 11, 18, 24, 25, 40, 45, 61, 62, 63, 72, 89, 108, 109, 110, 111, 112, 113]
    Apr = Float64[-6.748099080141577400e-001, 9.587587238674218900e-001, -3.380262965696901900e-001, -6.289006455522118100e-001, 1.458932368908371500e+000, 8.103505672584934100e-001, -8.308676110497204000e-002, 2.638152658681099300e+000, 1.000625209233638700e+000, 1.443285740856571400e-001, 1.198960671699122800e+000, 8.042866751960948700e-001, -2.146877760548437900e-002, 1.252109754253975600e+000, 1.681372772953815700e+000, 1.160558110039982400e+000, 8.971595941547876300e-001, 9.934462168535226300e-001, 4.030425643552826300e-003, 4.520907269135507900e-001, 1.957531681292445900e+000, 8.367513704018149100e-001, -8.919073879004054700e-001, 1.364968014685880000e-001, -5.391443069926015200e-001, 1.218605709809697200e+000, 1.379628849398701500e+000, 7.705722581573115600e-002, -9.721580354507602900e-001, 1.341091680254979900e+000, 1.416856131617781900e+000, 1.367782802004573000e+000, 9.310015370881947700e-001, 1.025960903569425500e+000, 4.418055920652105800e-001, 1.797608571927105100e-001, -1.566236421056189400e+000, -1.293790728779894400e-001, -2.548793507296202000e-001, -1.327158583802827700e+000, 8.378869252146687900e-001, 9.208088134425553900e-001, 8.954835779157467300e-001, 1.872488747783192400e+000, 1.607749082930002300e-001, 1.703097419135647200e+000, 1.485966470102646500e+000, 1.950727391193988900e+000, 1.261282074772785800e+000, 2.116080672704367100e+000, 4.418819171779880300e-001, -6.070414095655145600e-001, 2.236741804227308500e-001, -1.390565074056968600e+000, -3.689063294445240100e-001, 6.275925564016111600e-001, 9.419119351808997700e-001, 2.051844704916749900e-002, -2.055013180419911200e+000, -8.138411920140069500e-001, -1.670986440999406000e+000, -4.312826132928389000e-001, -1.317178055102502500e-001, 4.129877571668971700e-002, -1.178000809237306800e-001, -5.094809536848188300e-002, 4.887140456297478000e-001, -6.048920026011540600e-001, -1.622657687025348600e+000, 9.468924561461561400e-002, -1.076752006167445400e+000, -4.544402225817885700e-001, -1.138681599263542200e+000, 6.297099350763062400e-001, 7.161553122227841700e-001, -2.326213079050594200e-001, -2.421425730439526200e-001, -3.527400894253913700e-001, -5.991000491511071500e-001, -4.501982881436522000e-001, -1.073527618321097300e+000, -1.837831580649836900e+000, -1.093506922417250900e+000, 3.123080272148109900e-002, 4.764624397863112500e-001, -9.744593174286131400e-001, -5.475935551016516400e-001, -1.185577848308040100e+000, -1.512429413719079500e+000, -8.503358079552569800e-001, -1.048119449917124500e+000, -2.950353712891531400e-001, 4.436355670528358600e-001, -1.976143927829304500e+000, -1.154872945269262400e-001, -8.267783537066453200e-001, -2.524533500828647400e-001, -2.934013148063512300e-001, -9.559308197633347200e-001, -5.099545278978906400e-001, 7.104569782287112300e-001, 2.586151051423653000e-001, 9.670875274511733200e-001, -1.580472136196336100e-001, 4.036019169202506600e-001, -4.419329075671284400e-001, 1.335760493316310000e-001, -2.443367984060260100e-001, 5.346171843774407700e-001, 9.240837092450902800e-002, 1, -7.627629774659758300e+000, -8.317906940770274800e+000, -7.063111046330970300e+000, -8.991213209598194900e+000, -7.132056801314573300e+000, -6.933147988332642300e+000, -8.002850393982971200e+000, -7.763929473153173600e+000, -7.316196508551362500e+000, -6.857489440772187800e+000, -6.933435917031895100e+000, -8.453055934420948600e+000, -7.746956878849250200e+000, -6.987480243537588000e+000, -6.898256465718047200e+000, -7.035033952525910900e+000, -7.332300809407010800e+000, -8.868148466875196000e+000, -7.895445920224761800e+000, -8.477688226774738200e+000, -6.724790281248850500e+000, -8.026080227654171400e+000, -8.247796522362861500e+000, -7.893131787368933100e+000, -7.953077643849715300e+000, -7.085517647070380700e+000, -7.757918623860442500e+000, -9.404980440162075800e+000, -9.111021916635516900e+000, -6.753547752596220200e+000, -6.208328120779675900e+000, -7.460712402671073400e+000, -8.051154905293161700e+000, -6.979795864345947000e+000, -7.285373299414776400e+000, -7.935549231771511300e+000, -7.333663096710840700e+000, -7.549354579130950400e+000, -7.125439770244579300e+000, -8.444455837718113100e+000, -7.607705936847165300e+000, -6.947711222611268100e+000, -7.559561213742542300e+000, -6.633232534378238100e+000, -7.285067573640397100e+000, -7.198034090285221800e+000, -8.301173708858458500e+000, -7.740676955875905200e+000, -8.389317110487786200e+000, -7.030503781850090800e+000, 7.483837125022816400e-001, 9.972172897309477000e-001, 8.635272324715717000e-001, 3.437357445286000600e-001, 1.356057545748600200e+000, 7.747994541454063000e-001, -5.029337660688439900e-001, 7.320209587792132500e-001, 9.007328043370961300e-001, 1.839182838430090300e+000, -1.352604803382542600e+000, 1.440339792882617200e+000, 5.409160633286138600e-001, -9.096356667070149000e-002, 2.145906552675194100e-001, 1.749366326171780000e+000, 1.710862586429138200e+000, 5.617170924142114200e-001, 8.075937970441778900e-001, 5.294053117812258500e-001, 8.007328881544743000e-001, 6.824506852180559600e-001, -4.466422007073587200e-001, 4.125165857556476500e-001, 1.264599468598445900e+000, 1.195482116831558800e+000, 5.815701470999368200e-001, 1.480145121963046200e+000, 1.228974130124139900e+000, -9.887375252743877400e-001, 4.957232610676656800e-001, 2.201242074660017000e-001, -8.949779589135298300e-002, 4.463714750361750200e-001, 1.104458211701637000e+000, 1.146057798756675600e+000, 1.607064328082671100e+000, -3.342540671866551800e-001, -3.733869343080511400e-001, 1.376984209826379600e+000, 1.134842638638303900e+000, 1.578816025652644900e+000, 8.631161307174450800e-003, 6.628552272325566300e-001, -4.235551988694613600e-001, -5.931434773175525700e-002, 1.107744286910284700e+000, 1.308043758809049000e+000, 1.474220016986390400e+000, 6.094036504382410700e-001, 2.586151051423653000e-001, 7.690092056142631100e-001, -1.580472136196336100e-001, 2.410499789745967300e-001, -4.419329075671284400e-001, 6.725262868347948700e-001, -2.443367984060260100e-001, 9.310955439202097700e-001, 9.240837092450902800e-002, 9.816510048711224500e-001, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -5.0e-001, -5.0e-001, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1]
    b = Float64[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, -0.5]
    pwork = ECOS.ECOS_setup(
        n,
        m,
        p,
        l,
        ncones,
        q,
        0,
        Gpr,
        Gjc,
        Gir,
        Apr,
        Ajc,
        Air,
        c,
        h,
        b,
    )
    ECOS.unsafe_add_settings(pwork, Dict(:maxit => 1))
    @test ECOS.ECOS_solve(pwork) == ECOS.ECOS_MAXIT
    ECOS.ECOS_cleanup(pwork, 0)
    # Now test without max iterations
    pwork = ECOS.ECOS_setup(
        n,
        m,
        p,
        l,
        ncones,
        q,
        0,
        Gpr,
        Gjc,
        Gir,
        Apr,
        Ajc,
        Air,
        c,
        h,
        b,
    )
    @test ECOS.ECOS_solve(pwork) == ECOS.ECOS_OPTIMAL
    ECOS.ECOS_cleanup(pwork, 0)
    return
end

"""
Original problem:

Max   x
s.t.  x = 3
      8u + 9v <= 25
      -x + u + v = 0
      x free
      u <= 0
      v >= 0

ECOS form:

Min  -x
s.t.  x           == 3
     -x +  u +  v == 0
          8u + 9v <= 25
           u      <= 0
               -v <= 0
"""
function test_primal_infeasible()
    G = SparseArrays.sparse(Float64[0  8  9; 0  1  0; 0  0 -1])
    A = SparseArrays.sparse(Float64[1  0  0; -1  1  1])
    infeas_prob = ECOS.ECOS_setup(
        3,
        3,
        2,
        3,
        0,
        Int[],
        0,
        G.nzval,
        G.colptr .- 1,
        G.rowval .- 1,
        A.nzval,
        A.colptr .- 1,
        A.rowval .- 1,
        Float64[-1, 0, 0],
        Float64[25, 0, 0],
        Float64[3, 0],
    )
    @test ECOS.ECOS_solve(infeas_prob) == ECOS.ECOS_PINF
    ECOS.ECOS_cleanup(infeas_prob, 0)
    return
end

test_options()
test_primal_infeasible()


text:
module UnstaggeredMesh2D2P
export Mesh

struct Mesh
  x
  y
  p
  q
end

end


text:
# MATH OPERATIONS

for op in (:sin, :cos, :tan, :atan, :asin, :acos, :log, :log1p, :log10, :sqrt,
    :exp, :expm1, :erf, :erfc, :cbrt, :lgamma, :acosh, :cosh, :asinh, :sinh, :atanh, :tanh, :factorial)
    @eval Base.($(quot(op))){T}(x::AFAbstractArray{T}) = AFArray{T}(@cxx af::($op)(x.array))
end

Base.gamma{T}(x::AFAbstractArray{T}) = AFArray{T}(@cxx af::tgamma(x))

import Base: +, -, abs, .^

# TODO: add! using +=, etc.

.^{T,S}(a::AFAbstractArray{T}, b::AFAbstractArray{S}) = AFArray{af_promote(T,S)}(af_pow(a,b))
.^{T}(a::Real, b::AFAbstractArray{T}) = AFArray{T}(af_pow(a,b))
.^{T}(a::AFAbstractArray{T}, b::Real) = AFArray{T}(af_pow(a,b))

export root
root{T,S}(a::AFAbstractArray{T}, b::AFAbstractArray{S}) = AFArray{af_promote(T,S)}(af_root(a,b))
root{T}(a::Real, b::AFAbstractArray{T}) = AFArray{T}(af_root(a,b))
root{T}(a::AFAbstractArray{T}, b::Real) = AFArray{T}(af_root(a,b))

import Base: abs, min, max

#Abs
abs{T}(a::AFAbstractArray{T}) = AFArray{T}(af_abs(a))

#Max
function max(a::AFAbstractArray, b::AFAbstractArray)
    out = af_max(a,b)
    AFArray{backend_eltype(out)}(out)
end
function max(a::AFAbstractArray, b::Real)
    if typeof(b) <: Integer
        b = Float32(b)
    end
    out = af_max(a,b)
    AFArray{backend_eltype(out)}(out)
end
max(a::Real, b::AFAbstractArray) = max(b, a)

#Min
function min(a::AFAbstractArray, b::AFAbstractArray)
    out = af_min(a, b)
    AFArray{backend_eltype(out)}(out)
end
function min(a::AFAbstractArray, b::Real)
    if typeof(b) <: Integer
        b = Float32(b)
    end
    out = af_min(a,b)
    AFArray{backend_eltype(out)}(out)
end
min(a::Real, b::AFAbstractArray) = min(b, a)

#Negation
-{T}(a::AFSubArray{T}) = AFArray{T}(icxx"0-$a;")
-{T}(a::AFArray{T}) = AFArray{T}(af_neg(a))

#Logical ops
import Base: ==, .==, .>, .<, .>=, .<=, |, .!=

#Equals
.==(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = AFArray{Bool}(af_equals(a,b))

#Greater
.>(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = AFArray{Bool}(af_gt(a,b))
.>=(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = AFArray{Bool}(af_ge(a,b))

#Lesser
.<(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = AFArray{Bool}(af_lt(a,b))
.<=(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = AFArray{Bool}(af_le(a,b))

#Or
|(a::AFAbstractArray{Bool}, b::AFAbstractArray{Bool}) = AFArray{Bool}(af_bitor(a,b))

#not equals
.!=(a::AFAbstractArray, b::AFAbstractArray) = AFArray{Bool}(af_neq(a,b))

#Complex functions
import Base:complex, real, imag

complex{T}(a::AFAbstractArray{T}) = AFArray{Complex{T}}(af_complex(a))
complex{T}(a::AFAbstractArray{T}, b::AFAbstractArray{T}) = AFArray{Complex{T}}(af_complex(a,b))
real{T}(a::AFAbstractArray{Complex{T}}) = AFArray{T}(af_real(a))
imag{T}(a::AFAbstractArray{Complex{T}}) = AFArray{T}(af_imag(a))


text:
@testset "Test that we can use ForwardDiff" begin
    L0 = 1.0e-8
    KxStar = 1.0e-10
    f = 4
    Rtot = [100.0, 1000.0, 10.0, 10000.0]
    IgGC = [1.0, 0.1, 3.0]
    Kav = rand(3, 4) * 1.0e7
    Kav[1, 2] = 1.0e9

    @testset "randomBinding" begin
        func = x -> polyfc(L0, KxStar, f, x, IgGC, Kav).Rmulti
        out = ForwardDiff.gradient(func, Rtot)
        @test eltype(out) == Float64
        @test length(out) == length(Rtot)

        func = x -> polyfc(L0, KxStar, x, Rtot, IgGC, Kav).Rmulti
        out = ForwardDiff.derivative(func, Float64(f))
        @test typeof(out) == Float64
        @test out > 0.0

        func = x -> polyfc(x, KxStar, f, Rtot, IgGC, Kav).Rmulti
        out = ForwardDiff.derivative(func, L0)
        @test typeof(out) == Float64

        func = x -> polyfc(L0, x, f, Rtot, IgGC, Kav).Rmulti
        out = ForwardDiff.derivative(func, KxStar)
        @test typeof(out) == Float64
    end

    Cplx = [1 0 3; 2 2 0; 1 1 2; 4 0 0]
    Ctheta = [0.1, 0.3, 0.5, 0.1]
    Ltheta = L0 .* Ctheta

    @testset "complexBinding" begin
        for i = 1:4
            func = x -> sum(polyc(x, KxStar, Rtot, Cplx, Ctheta, Kav)[i])
            out = ForwardDiff.derivative(func, L0)
            @test typeof(out) == Float64

            func = x -> sum(polyc(L0, x, Rtot, Cplx, Ctheta, Kav)[i])
            out = ForwardDiff.derivative(func, KxStar)
            @test typeof(out) == Float64

            func = x -> sum(polyc(L0, KxStar, x, Cplx, Ctheta, Kav)[i])
            out = ForwardDiff.gradient(func, Rtot)
            @test eltype(out) == Float64
            @test length(out) == length(Rtot)

            func = x -> sum(polyc(L0, KxStar, Rtot, Cplx, x, Kav)[i])
            out = ForwardDiff.gradient(func, Ctheta)
            @test eltype(out) == Float64
            @test length(out) == length(Ctheta)

            func = x -> sum(polycm(KxStar, Rtot, Cplx, x, Kav)[i])
            out = ForwardDiff.gradient(func, Ltheta)
            @test eltype(out) == Float64
            @test length(out) == length(Ltheta)
        end
    end
end


text:
using LinearAlgebra
"""
        get_lagrange_basis(n = 1, varcoeff = false)

    Get Lagrange basis function of order `n`.
    """
function get_lagrange_basis(n = 1, varcoeff = false)
    # option varcoeff is not used actually
    global x
    global h
    E = Matrix(1I, n+1, n+1)
    symString = ""
    for i = 1:n
        symString *= ("a" * string(i - 1, base=10) * " ")
    end
    symString *= ("a" * string(n, base=10))

    L = SymPy.symbols(symString)

    basis = Array{SymPy.Sym}(undef, n+1)
    # p is a polynomial of degree n
    p = 0

    for i = 1:n+1
        p = p + L[i]*x^(i-1)
    end
    points = []
    for i = 0:n
        push!(points, i * h / n)
    end

    for i = 1:(n + 1)
        Leq = zeros(Sym, n+1)

        for j = 1:(n + 1)
            eq = SymPy.subs(p, x, points[j])
            Leq[j] = eq
        end

        Leq = Leq .- E[1:(n+1), i]
        if n == 0
            L2 = zeros(Sym, 1)
            L2[1] = L
            L = L2
        end
        LL = SymPy.solve(Leq, collect(L))
        LL = convert(Dict{SymPy.Sym, SymPy.Sym}, LL)
        q = 0
        for j = 1:(n+1)
            q = q + LL[L[j]] * x^(j - 1)
        end
        basis[i] = q
    end
    basis
end

"""
    get_lagrange_em(p = 1, m = 0, n = 0)

Get Lagrange finite elements elementary matrices.

# Arguments
  * `p`: degree of polynomials in the basis.
  * `m`: number of derivatives on the first function.
  * `n`: number of derivatives on the second function.
"""
function get_lagrange_em(p = 1, m = 0, n = 0)
    global x
    global h
    M = Array{SymPy.Sym}(undef, p+1, p+1)
    F = get_lagrange_basis(p)
    for i = 1:p+1
        for j = 1:p+1
            M[i, j] = simplify(integrate(diff(F[i], x, m) * diff(F[j], x, n), (x, 0, h)))
        end
    end
    M
end

"""
    get_lagrange_em_varcoeff(p = 1, m = 0, n = 0, f = 1)

Get Hermite finite elements elementary matrices for variable coefficients.

# Arguments
* `p`: degree of polynomials in the basis.
* `m`: number of derivatives on the first function.
* `n`: number of derivatives on the second function.
* `f`: the variable coefficient.
"""
function get_lagrange_em_varcoeff(p = 1, m = 0, n = 0, f = 1)
    # does not work properly
    # TODO: fix it
    global x
    global h
    global xa
    global xb
    M = Array{SymPy.Sym}(undef, p+1, p+1)
    F = get_lagrange_basis(p, true)
    for i = 1:p+1
        for j = 1:p+1
            M[i, j] = simplify(integrate(diff(F[i], x, m) * diff(F[j], x, n) * f, (x, xa, xb)))
            M[i, j] = simplify(subs(M[i, j], h, xb - xa))
        end
    end
    M
end

"""
    get_square_lagrange_em((px, py) = (1, 1), (mx, my) = (0, 0), (nx, ny) = (0, 0))

Get Lagrange finite elements elementary matrices for a squared element.

# Arguments
  * `(px, py)` : degree of polynomials in the basis.
  * `(mx, my)` : number of derivatives on the first function wrt x and y
  * `(nx, ny)` : number of derivatives on the second function wrt x and y

# Remarks
  The current version works only for px, py in {1, 2}
"""
function get_square_lagrange_em((px, py) = (1, 1), (mx, my) = (0, 0), (nx, ny) = (0, 0))
    Mx = get_lagrange_em(px, mx, nx)
    My = get_lagrange_em(py, my, ny)

    node_x = [(1, 2, 2, 1),
              (1, 3, 3, 1, 2, 3, 2, 1, 2)]
    node_y = [(1, 1, 2, 2),
              (1, 1, 3, 3, 1, 2, 3, 2, 2)]

    p = (px + 1)*(py + 1)
    M = Array{SymPy.Sym}(undef, p, p)

    for i = 1:p
        for j = 1:p
            M[i,j] = Mx[node_x[px][i], node_x[px][j]] * My[node_y[py][i], node_y[py][j]]
        end
    end
    M
end


text:
"""
    Unexported function to deal with the coordinates required by VoronoiDelaunay
"""
function _nrm(v, m, M, n, N)
    vn = (v .- m) ./ (M - m)
    return (vn .* (N - n)) .+ n
end

"""
    wombling(x::Vector{T}, y::Vector{T}, z::Vector{T}) where {T<:Number}

Wrapper function that implements the triangulation wombling algorithm for points
that are irregularly arranged in space.
"""
function wombling(x::Vector{T}, y::Vector{T}, z::Vector{T}) where {T<:Number}
    length(x) >= 3 || throw(DimensionMismatch("x must have a minimum length of 4"))
    length(x) == length(y) ||
        throw(DimensionMismatch("x and y must have the same dimension"))
    length(x) == length(z) ||
        throw(DimensionMismatch("x and z must have the same dimension"))

    # Get the range of values for the points
    min_value, max_value = extrema(vcat(x, y))

    # Project the points in the correct range of values
    nx = _nrm(x, min_value, max_value, min_coord, max_coord)
    ny = _nrm(y, min_value, max_value, min_coord, max_coord)
    px = Point2D[VoronoiDelaunay.Point(nx[i], ny[i]) for i in eachindex(nx)]

    # Build the Delaunay triangulation
    tess = DelaunayTessellation()
    sizehint!(tess, length(px))

    # This bit is important because we will need to find the correct index of each of these things
    for p in px
        push!(tess, p)
    end
    triangles = unique(tess)

    _M = zeros(T, length(triangles))
    _θ = similar(_M)
    _X = similar(_M)
    _Y = similar(_M)

    # Get the rate of change for the points
    for (i, triangle) in enumerate(triangles)
        centroid_x = mean([triangle._a._x, triangle._b._x, triangle._c._x])
        centroid_y = mean([triangle._a._y, triangle._b._y, triangle._c._y])
        simplices_coordinates = [
            findfirst(p -> (p._x == triangle._a._x)&(p._y == triangle._a._y), px),
            findfirst(p -> (p._x == triangle._b._x)&(p._y == triangle._b._y), px),
            findfirst(p -> (p._x == triangle._c._x)&(p._y == triangle._c._y), px)
        ]
        _x = x[simplices_coordinates]
        _y = y[simplices_coordinates]
        _z = z[simplices_coordinates]
        _M[i], _θ[i] = SpatialBoundaries._rateofchange(_x, _y, _z)
        _X[i] = mean(_x)
        _Y[i] = mean(_y)
    end

    # Rate of change and direction
    return TriangulationWomble(_M, _θ, _X, _Y)
end

"""
    wombling(x::Vector{T}, y::Vector{T}, z::Matrix{T}) where {T<:Number}

Wrapper function that implements the lattice wombling algorithm for points that
are regularly arranged in space. Note that the matrix is presented in a way that
is flipped, *i.e.* the `x` coordinates corresponds to the rows, and the `y`
coordinates correspond to the columns. If you want to think of `x` and `y` as
geographic coordinates, `y` are the longitudes, and `x` are the latitudes. Using
the bindings for `SimpleSDMLayers`, this conversion will be performed
automatically.
"""
function wombling(x::Vector{T}, y::Vector{T}, z::Matrix{T}) where {T<:Number}
    issorted(x) || throw(ArgumentError("The values of x must be sorted and increasing"))
    issorted(y) || throw(ArgumentError("The values of y must be sorted and increasing"))
    length(x) == size(z, 1) || throw(
        DimensionMismatch("The length of x must be equal to the first dimension of z")
    )
    length(y) == size(z, 2) || throw(
        DimensionMismatch("The length of y must be equal to the second dimension of z")
    )

    _M = zeros(T, size(z) .- 1)
    _θ = similar(_M)

    for j in 1:size(_M, 2), i in 1:size(_M, 1) # womble along a 2x2 window
        tmp = z[i:(i + 1), j:(j + 1)]
        _M[i, j], _θ[i, j] = SpatialBoundaries._rateofchange(
            x[i:(i + 1)], y[j:(j + 1)], tmp
        )
    end

    # Rate of change and direction
    return LatticeWomble(
        _M,
        _θ,
        x[1:(end - 1)] .+ 0.5 .* vec(x[2:end] .- x[1:(end - 1)]),
        y[1:(end - 1)] .+ 0.5 .* vec(y[2:end] .- y[1:(end - 1)]),
    )
end

"""
    wombling(m::Matrix{T}) where {T<:Number}

Shortcut to womble a matrix (using lattice wombling) when no x and y positions
are given - the cell size in each dimension is expected to be 1.
"""
function wombling(m::Matrix{T}) where {T<:Number}
    return wombling(convert.(T, 1:size(m, 1)), convert.(T, 1:size(m, 2)), m)
end


text:
# apollonian gasket constructed in sketch environment
# please dont` use large subdivision levels if you cant` afford supercomputer
#

using ShapeFactory

app = dispatch("CATIA.Application")
prt = app.ActiveDocument.Part

xy = prt.FindObjectByName("xy plane")

gst = prt.HybridBodies.Add()
gst.Name = "Output"

sk1 = gst.HybridSketches.Add(xy)
f2d = sk1.OpenEdition()
con = sk1.Constraints

function base_circle(rad)
    pt0 = f2d.CreatePoint(0, 0)
    con.AddMonoEltCst.(0, pt0)

    ax0 = f2d.CreateLine(0, 0, 0, rad * 3)
    ax0.StartPoint = pt0
    ax0.Construction = true
    con.AddMonoEltCst(13, ax0)

    c0 = f2d.CreateClosedCircle(0, 0, rad)
    c0.CenterPoint = pt0
    rad = con.AddMonoEltCst.(14, c0)
    rad.Dimension.ValuateFromString("$(rad)mm")

    return c0, ax0
end

function base_triad(cc0, ax0)

    cx = f2d.CreateClosedCircle(0, 0, cc0.Radius - 1)
    cx.CenterPoint = cc0.CenterPoint
    cx.Construction = true

    tri = map((2pi/3):(2pi/3):2pi) do ang
        pt = f2d.CreatePoint(sin(ang) * cx.Radius, cos(ang)*cx.Radius)
        pt.Construction = true
        cc = f2d.CreateClosedCircle(sin(ang), cos(ang), 1)
        cc.CenterPoint = pt
        rad = con.AddMonoEltCst.(14, cc)
        rad.Dimension.ValuateFromString("$(cc0.Radius / 2)mm")
        cc => rad
    end

    tri, rad = first.(tri), last.(tri)

    con.AddBiEltCst.(2, getproperty.(tri, :CenterPoint), cx)

    push!(tri, tri[1])
    for n in 1:3
        con.AddBiEltCst(4, tri[n], tri[n+1])
    end

    con.AddBiEltCst(2, tri[1].CenterPoint, ax0)
    map(x -> x.Deactivate(), rad)
    con.AddBiEltCst.(4, tri[1:end-1], cc0)

    return tri[1:end-1]
end

function circle_tri(a, b, c)
    loop = (a, b, c, a)
    cntr = [0, 0, 0]
    for n in 1:3
        pt    = f2d.CreatePoint(0, 0)
        pc    = con.AddBiEltCst.(2, pt, (loop[n], loop[n+1]))
        cntr += pt.GetCoordinates([0, 0, 0])[2] ./ 3
        map(x -> x.Deactivate(), pc)
    end
    cir = f2d.CreateClosedCircle(cntr[1], cntr[2], 0.01)
    con.AddBiEltCst.(4, cir, sort([a, b, c], by=x -> x.Radius, rev=true))
    return cir
end

function subdivide(a, b, c, stop)
    loop = (a, b, c, a)
    x = circle_tri(a, b, c)
    if stop > 1
        for n in 1:3
            subdivide(loop[n], loop[n+1], x, stop-1)
        end
    end
end

c0, ax = base_circle(10)
a, b, c = base_triad(c0, ax)
subdivide(a, b, c, 2)

loop = (a, b, c, a)
for n in 1:3
    subdivide(loop[n], loop[n+1], c0, 3)
end

sk1.CloseEdition()
prt.Update()

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory


text:
"""
    docstring_ext

Internal function to help generate the docstring for a package. See
[`metadata_model`](@ref).
"""
function docstring_ext(T; descr::String="")
    package_name = MLJModelInterface.package_name(T)
    package_url  = MLJModelInterface.package_url(T)
    model_name   = MLJModelInterface.name(T)
    # the message to return
    message      = "$descr"
    message     *= "\n→ based on [$package_name]($package_url)."
    message     *= "\n→ do `@load $model_name pkg=\"$package_name\"` to " *
                   "use the model."
    message     *= "\n→ do `?$model_name` for documentation."
end

"""
    metadata_pkg(T; args...)

Helper function to write the metadata for a package providing model `T`.
Use it with broadcasting to define the metadata of the package providing
a series of models.

## Keywords

* `name="unknown"`   : package name
* `uuid="unknown"`   : package uuid
* `url="unknown"`    : package url
* `julia=missing`    : whether the package is pure julia
* `license="unknown"`: package license
* `is_wrapper=false` : whether the package is a wrapper

## Example

```
metadata_pkg.((KNNRegressor, KNNClassifier),
    name="NearestNeighbors",
    uuid="b8a86587-4115-5ab1-83bc-aa920d37bbce",
    url="https://github.com/KristofferC/NearestNeighbors.jl",
    julia=true,
    license="MIT",
    is_wrapper=false)
```
"""
function metadata_pkg(T; name::String="unknown",
                         uuid::String="unknown",
                         url::String="unknown",
                         julia::Union{Missing,Bool}=missing,
                         license::String="unknown",
                         is_wrapper::Bool=false)
    ex = quote
        MLJModelInterface.package_name(::Type{<:$T})    = $name
        MLJModelInterface.package_uuid(::Type{<:$T})    = $uuid
        MLJModelInterface.package_url(::Type{<:$T})     = $url
        MLJModelInterface.is_pure_julia(::Type{<:$T})   = $julia
        MLJModelInterface.package_license(::Type{<:$T}) = $license
        MLJModelInterface.is_wrapper(::Type{<:$T})      = $is_wrapper
    end
    parentmodule(T).eval(ex)
end

"""
    metadata_model(`T`; args...)

Helper function to write the metadata for a model `T`.

## Keywords

* `input=Unknown` : allowed scientific type of the input data
* `target=Unknown`: allowed sc. type of the target (supervised)
* `output=Unknown`: allowed sc. type of the transformed data (unsupervised)
* `weights=false` : whether the model supports sample weights
* `descr=""`      : short description of the model
* `path=""`       : where the model is (usually `PackageName.ModelName`)

## Example

```
metadata_model(KNNRegressor,
    input=MLJModelInterface.Table(MLJModelInterface.Continuous),
    target=AbstractVector{MLJModelInterface.Continuous},
    weights=true,
    descr="K-Nearest Neighbors classifier: ...",
    path="NearestNeighbors.KNNRegressor")
```
"""
function metadata_model(T; input=Unknown,
                           target=Unknown,
                           output=Unknown,
                           weights::Bool=false,
                           descr::String="",
                           path::String="")
    if isempty(path)
        pname = MLJModelInterface.package_name(T)
        mname = MLJModelInterface.name(T)
        path = "MLJModels.$(pname)_.$(mname)"
    end
    ex = quote
        MLJModelInterface.input_scitype(::Type{<:$T})    = $input
        MLJModelInterface.output_scitype(::Type{<:$T})   = $output
        MLJModelInterface.target_scitype(::Type{<:$T})   = $target
        MLJModelInterface.supports_weights(::Type{<:$T}) = $weights
        MLJModelInterface.load_path(::Type{<:$T})        = $path

        MLJModelInterface.docstring(::Type{<:$T}) =
            MLJModelInterface.docstring_ext($T; descr=$descr)
    end
    parentmodule(T).eval(ex)
end


text:
# This script actually solves the RBCCampbellCochrane model with a risk-adjusted linearization
# and times the methods, if desired
using BenchmarkTools, RiskAdjustedLinearizations, Test
include("rbc_cc.jl")

# Settings: what do you want to do?
time_methods        = false
numerical_algorithm = :relaxation
autodiff            = false
n_strips            = 3

# Set up
autodiff_method = autodiff ? :forward : :central
m_rbc_cc = RBCCampbellCochraneHabits()
m = rbc_cc(m_rbc_cc, n_strips)
z0 = copy(m.z)
y0 = copy(m.y)
Ψ0 = copy(m.Ψ)

# Solve!
solve!(m; algorithm = numerical_algorithm, autodiff = autodiff_method)

if n_strips == 0
    if numerical_algorithm == :relaxation
        sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "rbccc_sss_iterative_output.jld2"), "r")
    elseif numerical_algorithm == :homotopy
        sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "rbccc_sss_homotopy_output.jld2"), "r")
    end
elseif n_strips == 3
    if numerical_algorithm == :relaxation
        sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "rbccc_sss_iterative_N3_output.jld2"), "r")
    elseif numerical_algorithm == :homotopy
        sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "rbccc_sss_homotopy_N3_output.jld2"), "r")
    end
end

if n_strips in [0, 3]
    @test isapprox(sssout["z_rss"], m.z, atol=1e-4)
    @test isapprox(sssout["y_rss"], m.y, atol=2e-4)
    @test isapprox(sssout["Psi_rss"], m.Ψ, atol=1e-4)
end

if time_methods
    println("Deterministic steady state")
    @btime begin
        solve!(m, z0, y0; algorithm = :deterministic, autodiff = autodiff_method, verbose = :none)
    end

    # Use deterministic steady state as guesses
    solve!(m, z0, y0; algorithm = :deterministic, autodiff = autodiff_method, verbose = :none)
    zdet = copy(m.z)
    ydet = copy(m.y)
    Ψdet = copy(m.Ψ)

    println("Relaxation method")
    @btime begin # called the "iterative" method in the original paper
        solve!(m, zdet, ydet, Ψdet; algorithm = :relaxation, autodiff = autodiff_method, verbose = :none)
    end

    println("Relaxation method with automatic differentiation")
    @btime begin # called the "iterative" method in the original paper
        solve!(m, zdet, ydet, Ψdet; algorithm = :relaxation, autodiff = :forward, verbose = :none)
    end

    println("Relaxation method with Anderson acceleration")
    @btime begin # called the "iterative" method in the original paper
        solve!(m, zdet, ydet, Ψdet; algorithm = :relaxation, use_anderson = true, m = 3, autodiff = autodiff_method, verbose = :none)
    end

    println("Homotopy method")
    @btime begin # called the "continuation" method in the original paper, but is called homotopy in the original code
        solve!(m, zdet, ydet, Ψdet; algorithm = :homotopy, autodiff = autodiff_method, verbose = :none)
    end
end


text:
using Test

const _TUTORIAL_DIR = joinpath(@__DIR__, "..", "docs", "src", "tutorials")

@testset "run_examples.jl" begin
    @testset "$dir" for dir in readdir(_TUTORIAL_DIR)
        for filename in readdir(joinpath(_TUTORIAL_DIR, dir))
            if filename[1] == '_'
                continue
            end
            path = joinpath(_TUTORIAL_DIR, dir, filename)
            include(path)
        end
    end
    @testset "Chordal" begin
        include("chordal_sparsity.jl")
        include("chordal_sparsity_with_domain.jl")
    end
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Test, RubiSymbolics
using RubiSymbolics: @variables, simplify, isequal, expand, Elliptic, HypergeometricFunctions, @rule, Fixpoint, Postwalk, Chain

const print_incorrect_tests = false

r1 = @rule (~base::(x -> x isa Number && x<0))^(~exp::(x -> !(x isa Integer))) --> Complex(~base)^~exp
r1 = @rule (-(~base::(x -> x isa Number && x>0)))^(~exp::(x -> !(x isa Integer))) --> Complex(-(~base))^~exp
exp_fixer = Fixpoint(Postwalk(Chain([r1])))
function preprocess(expr::Expr)
    return exp_fixer(expr)
end

macro test_int(expr)
    if (expr.head == :vect && length(expr.args) == 4)
        preprocessed = preprocess(expr)
        (integrand, symbol, steps, result) = preprocessed.args

        # return :(@test isequal(simplify(integrate($(esc(integrand)), $(esc(symbol)))), simplify($(esc(result)))))
        return :(@test isequal(integrate($(esc(integrand)), $(esc(symbol))), $(esc(result))))
        # quote
        #     isequal(simplify(integrate($(esc(integrand)), $(esc(symbol)))), simplify($(esc(result))))
        #     @test true
        # end
    else
        if print_incorrect_tests
            println("Incorrect test problem at ", __source__)
        end
        return :(@test_broken false)
    end
end

# include(".jl")
include("/home/main/.julia/dev/RubiSymbolics/test/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/.jl")


text:
module MeshIO

using GeometryBasics
using ColorTypes
using Printf
import FileIO

using GeometryBasics: raw, value, decompose_normals, convert_simplex
import FileIO: DataFormat, @format_str, Stream, File, filename, stream
import FileIO: skipmagic, add_format

import Base.show

include("io/off.jl")
include("io/ply.jl")
include("io/stl.jl")
include("io/obj.jl")
include("io/2dm.jl")
include("io/msh.jl")

"""
    load(fn::File{MeshFormat}; pointtype=Point3f0, uvtype=Vec2f0,
         facetype=GLTriangleFace, normaltype=Vec3f0)

"""
function load(fn::File{format}; element_types...) where {format}
    open(fn) do s
        skipmagic(s)
        load(s; element_types...)
    end
end

function save(fn::File{format}, msh::AbstractMesh) where {format}
    open(fn, "w") do s
        save(s, msh)
    end
end

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end # module


text:
using StructuresKit

E = 29000.0
σy = 50.0                   #steel yield stress
σy1 = 50.0                  #steel yield stress
σu = 65.3                   #steel ultimate stress
σf=σy * (1 + 0.10)      #steel fracture stress#
ϵy = σy / E              #steel yield strain
ϵy1 = σy / E * 10    #steel strain at end of yield plateau
ϵu = 0.18             #steel ultimate strain
ϵf = 0.21             #steel fracture strain
n = 11

σ, ϵ  = MaterialModels.steel(σy, σy1, σu, σf, ϵy, ϵy1, ϵu, ϵf, n)

using Plots
plot(σ, ϵ, markershape = :o)

shape_name = "W14X90"

shape_info = CrossSection.AISC(shape_name)

n_Wshape=(4, 2, 4, 4, 4)

xcoords, ycoords = CrossSection.wshape_nodes(shape_info, n_Wshape)

plot(xcoords, ycoords, markershape = :o)

# using TriangleMesh

mesh_size = 0.01
mesh = CrossSection.triangular_mesh(xcoords, ycoords, mesh_size)

Ai, cxi, cyi = CrossSection.triangular_mesh_properties(mesh)

cx, cy = CrossSection.centroid_from_cells(Ai, cxi, cyi)

Ix = CrossSection.moment_of_inertia_from_cells(Ai, cyi, cy)
Iy = CrossSection.moment_of_inertia_from_cells(Ai, cxi, cx)

# #discretize cross-section with a triangular mesh
# function triangular_mesh(xcoords, ycoords, mesh_size)

#     num_nodes = length(xcoords)
#     num_segments = num_nodes

#     # n_point, n_point_marker, n_point_attribute, n_segment, n_holes
#     poly = TriangleMesh.Polygon_pslg(num_nodes, 1, 0, num_segments, 0)

#     node = [xcoords ycoords]
#     set_polygon_point!(poly, node)

#     node_marker = ones(Int, num_nodes, 1)
#     set_polygon_point_marker!(poly, node_marker)

#     segments = zeros(Int, num_segments, 2)
#     for i=1:num_segments

#         if i == num_segments
#             segments[i, 1:2] = [i, 1]
#         else
#             segments[i, 1:2] = [i, i+1]
#         end

#     end

#     set_polygon_segment!(poly, segments)

#     segment_markers = ones(Int, num_segments)
#     set_polygon_segment_marker!(poly, segment_markers)

#     #switches from https://www.cs.cmu.edu/~quake/triangle.html
#     switches = "penvVa" * string(mesh_size) * "D"

#     mesh = create_mesh(poly, switches)

#     return mesh

# end

#calculate triangle areas

# # https://keisan.casio.com/has10/SpecExec.cgi?path=05000000.Mathematics%252F01000500.Plane%2520geometry%252F10010300.Area%2520of%2520a%2520triangle%2520with%2520three%2520points%252Fdefault.xml&charset=utf-8
# function triangle_area(x1, y1, x2, y2, x3, y3)

#     A = abs((x1*y2 + x2*y3 + x3*y1 - y1*x2 - y2*x3 - y3*x1)/2)

# end

# #https://www.mathopenref.com/coordcentroid.html
# function triangle_centroid(x1, y1, x2, y2, x3, y3)

#     cx = (x1 + x2 + x3)/3
#     cy = (y1 + y2 + y3)/3

#     return cx, cy

# end

# function triangulation_properties(mesh)

#     #calculate cell area and centroid
#     Ai = zeros(Float64, mesh.n_cell)
#     cxi = zeros(Float64, mesh.n_cell)
#     cyi = zeros(Float64, mesh.n_cell)

#     for i = 1:mesh.n_cell

#         p1 = mesh.cell[1, i]
#         p2 = mesh.cell[2, i]
#         p3 = mesh.cell[3, i]

#         x1 = mesh.point[1, p1]
#         y1 = mesh.point[2,p1]
#         x2 = mesh.point[1, p2]
#         y2 = mesh.point[2,p2]
#         x3 = mesh.point[1, p3]
#         y3 = mesh.point[2,p3]

#     Ai[i] = triangle_area(x1, y1, x2, y2, x3, y3)

#     cxi[i], cyi[i] = triangle_centroid(x1, y1, x2, y2, x3, y3)

#     end

#     return Ai, cxi, cyi

# end

# #calculate cross-sectional area from cells

# # function area_from_cells(Ai)

# #     A = sum(Ai)

# # end

# # #calculate cross-section centroid from cells

# # function centroid_from_cells(Ai, cxi, cyi)

# #     A = area_from_cells(Ai)

# #     cx = sum(cxi .* Ai) / A
# #     cy = sum(cyi .* Ai) / A

# #     return cx, cy

# # end

# cx, cy = CrossSection.centroid_from_cells(Ai, cxi, cyi)

# # function moment_of_inertia_from_cells(Ai, ci, c)

# #     I = sum(((c .- ci) .^2 .* Ai))

# #     return I

# # end

# Ix = moment_of_inertia_from_cells(Ai, cyi, cy)
# Iy = moment_of_inertia_from_cells(Ai, cxi, cx)

# #calculate triangle centroids

# # for i = 1:mesh.n_cell

# #     p1 = mesh.cell[1, i]
# #     p2 = mesh.cell[2, i]
# #     p3 = mesh.cell[3, i]

# #     x1 = mesh.point[1, p1]
# #     y1 = mesh.point[2,p1]
# #     x2 = mesh.point[1, p2]
# #     y2 = mesh.point[2,p2]
# #     x3 = mesh.point[1, p3]
# #     y3 = mesh.point[2,p3]

# #    cx[i] = triangle_area(x1, y1, x2, y2, x3, y3)

# # end

# zcoords = zeros(Float64, mesh.n_point)
# coordinates = [mesh.point[1,:] mesh.point[2,:] zcoords]

# connectivity = mesh.cell'

# using Makie
# scene = Makie.poly(coordinates, connectivity, color=:lightgray, shading=:true, show_axis=:true, overdraw=:false, strokecolor = (:black, 0.6), strokewidth = 1)

# poly = polygon_Lshape()
# mesh = create_mesh(poly, info_str="my mesh", voronoi=true, delaunay=true, set_area_max=true)

# node = [1.0 0.0 ; 0.0 1.0 ; -1.0 0.0 ; 0.0 -1.0 ;
#         0.25 0.25 ; -0.25 0.25 ; -0.25 -0.25 ; 0.25 -0.25]

#         # size is number_segments x 2
# seg = [1 2 ; 2 3 ; 3 4 ; 4 1 ; 5 6 ; 6 7 ; 7 8 ; 8 5]

# # all points get marker 1
# node_marker = [ones(Int,4,1) ; 2*ones(Int,4,1)]
# # last segment gets a different marker
# seg_marker = [ones(Int,4) ; 2*ones(Int,4)]

# # size is number_points x number_attr
# node_attr = rand(8,2)

# # size is number_holes x 2
# hole = [0.5 0.5]

# polygon_struct_from_points

# shape_mesh, cross_section_edges = Mesh.open_cross_section_tessellation(xcoords, ycoords)

# connectivity = zeros(Int, length(shape_mesh), 3)
# for i=1:length(shape_mesh)
#     connectivity[i, :] = shape_mesh[i]
# end


text:
using PList

function completion_choices_plist(
    filepath::String,
    row::Int,
    col::Int,
    gocode::String)

    document = open(readlines, filepath)
    # byte offset of cursor position from the beginning of file
    cursor = sum(map(length, document[1:row - 1])) + col - 1
    (stream, process) = readsfrom(`$gocode -f=csv -in=$filepath autocomplete $cursor`)
    suggestion_table = readcsv(stream)
    no_rows = size(suggestion_table)[1]
    choices = map(make_completion_choices, {suggestion_table[row,:] for row in 1:no_rows})
    writeplist_string(choices)
    # suggestions = writeplist_string({{"display" => row} for row in suggestion_table[:,3]})
end

function make_completion_choices(comp)
    image = comp[1]
    match = comp[3]
    args = comp[end]
    insert = snippet_from_func_args(args)
    display = string(match, image == "func" ? replace(args, "func", "") : " " * args)

    { "match" => match, "display" =>  display, "insert" => insert, "image" => image }
end


text:
# Implementation of functionals (i.e. higher-order functions).

# Implementation of sensitivities w.r.t. `map`.

# Build `@explicit_intercepts`-like calls for `map` with a variable number of arguments.
# We set an arbitrary cutoff of 10 input arrays, which should be sufficient, as mapping
# over more than 3 arrays simultaneously is (anecdotally) pretty uncommon. Imposing this
# cutoff alleviates the need to overwrite a Base method to call `invoke` due to method
# ambiguities, though even that can be insufficient to resolve them, as in issue #136.
const ArrayOrNode = Union{AbstractArray, Node{<:AbstractArray}}
for nargs = 1:10
    name = Symbol(:A, nargs)
    args = Expr[:($(Symbol(:A, i))::$ArrayOrNode) for i = 1:nargs-1]
    push!(args, :($name::Node{<:AbstractArray}))
    @eval function Base.map(f, $(args...), As::$ArrayOrNode...)
        Branch(map, (f, $(map(i->Symbol(:A, i), 1:nargs)...), As...), getfield($name, :tape))
    end
end

# Compute sensitivity w.r.t. the N^{th} input, N > 1.
∇(::typeof(map), ::Type{Arg{N}}, p, y, ȳ, f::Function, A::∇Array...) where N =
    _∇(map, Arg{N-1}, p, y, ȳ, f, A...)
_∇(::typeof(map), arg::Type{Arg{N}}, p, y, ȳ, f::Function, A::∇Array...) where N =
        map((ȳn, An...)->ȳn * fmad(f, An, Val{N}), ȳ, A...)

# Implementation of sensitivities w.r.t. `broadcast`.
using Base.Broadcast
using Base.Broadcast: Broadcasted, broadcastable, broadcast_axes, combine_axes, result_style

struct NodeStyle{S} <: BroadcastStyle end

Base.BroadcastStyle(::Type{<:Node{T}}) where {T} = NodeStyle{BroadcastStyle(T)}()

Base.BroadcastStyle(::NodeStyle{S}, ::NodeStyle{S}) where {S} = NodeStyle{S}()
function Base.BroadcastStyle(::NodeStyle{S1}, ::NodeStyle{S2}) where {S1,S2}
    promoted = result_style(S1, S2)
    promoted isa Broadcast.Unknown ? promoted : NodeStyle{promoted}()
end
function Base.BroadcastStyle(::NodeStyle{S}, B::BroadcastStyle) where {S}
    promoted = result_style(S, B)
    promoted isa Broadcast.Unknown ? promoted : NodeStyle{promoted}()
end

Base.axes(x::Node) = axes(unbox(x))
Broadcast.broadcastable(x::Node) = x

# eagerly construct a Branch when encountering a Node in broadcasting
function Broadcast.broadcasted(::NodeStyle, f, args...)
    tape = getfield(args[findfirst(x -> x isa Node, args)], :tape)
    return Branch(broadcast, (f, args...), tape)
end

"""
    broadcastsum!(f::Function, add::Bool, z, As...)

Broadcast f over As and reduce to z by summing. If add is true, then the result is added to
the current value of z, otherwise it is overwritten.
"""
function broadcastsum!(f::Function, add::Bool, z, As...)
    tmp_shape = map(length, combine_axes(As...))
    if size(z) != tmp_shape
        tmp = Array{eltype(z)}(undef, tmp_shape)
        return sum!(z, broadcast!(f, tmp, As...), init=!add)
    else
        return add ?
            broadcast!((z, x...)->z + f(x...), z, z, As...) :
            broadcast!(f, z, As...)
    end
end

"""
    broadcastsum(f, add::Bool, z::AbstractArray, As...)

Allocating version of broadcastsum! specialised for Arrays.
"""
broadcastsum(f, add::Bool, z::AbstractArray, As...) =
    broadcastsum!(f, add, Array{eltype(z)}(undef, size(z)), As...)

"""
    broadcastsum(f, add::Bool, z::Number, As...)

Specialisation of broadcastsum to Number-sized outputs.
"""
function broadcastsum(f, add::Bool, z::Number, As...)
    tmp = Array{eltype(z)}(undef, map(length, combine_axes(As...)))
    return sum(broadcast!(f, tmp, As...)) + (add ? z : zero(z))
end

broadcastsum(f, add::Bool, z::Ref{<:Number}, As...) = broadcastsum(f, add, z[], As...)

# Compute sensitivity w.r.t. the N^{th} input, N > 1.
∇(::typeof(broadcast), ::Type{Arg{N}}, p, y, ȳ, f, A...) where N =
    _∇(broadcast, Arg{N-1}, p, y, ȳ, f, A...)
_∇(::typeof(broadcast), ::Type{Arg{N}}, p, y, ȳ, f, A...) where N =
        broadcastsum((ȳn, xn...)->ȳn * fmad(f, xn, Val{N}), false, A[N], ȳ, A...)

# Division from the right by a scalar.
import Base: /
@eval @explicit_intercepts $(Symbol("/")) Tuple{∇Array, ∇Scalar}
@inline ∇(::typeof(/), ::Type{Arg{1}}, p, z, z̄, x::∇Scalar, y::∇Array) =
    ∇(broadcast, Arg{2}, p, z, z̄, /, x, y)
@inline ∇(::typeof(/), ::Type{Arg{2}}, p, z, z̄, x::∇Scalar, y::∇Array) =
    ∇(broadcast, Arg{3}, p, z, z̄, /, x, y)

# Division from the left by a scalar.
import Base: \
@eval @explicit_intercepts $(Symbol("\\")) Tuple{∇Scalar, ∇Array}
@inline ∇(::typeof(\), ::Type{Arg{1}}, p, z, z̄, x::∇Array, y::∇Scalar) =
    ∇(broadcast, Arg{2}, p, z, z̄, \, x, y)
@inline ∇(::typeof(\), ::Type{Arg{2}}, p, z, z̄, x::∇Array, y::∇Scalar) =
    ∇(broadcast, Arg{3}, p, z, z̄, \, x, y)


text:
################################################################################
#
#   PLOTTING SITES
#
################################################################################

# PLOTTING SITES IN 2D (but complex)
function plotSiteComplex(
            site    :: S,
            radius  :: Real,
            color   :: Vector{<:Integer}
            ;
            detail  :: Int64 = 21,
            kwargs...
        ) where {L,S<:AbstractSite{L,2}}

    # get the coordinates
    x = point(site)[1]
    y = point(site)[2]
    # Make data of sphere
    n = detail
    u = range(0,stop=2*pi,length=n);
    px = cos.(u);
    py = sin.(u);
    px = (px.*radius) .+ x
    py = (py.*radius) .+ y
    # Plot the surface of the sphere
    PyPlot.fill(px,py, color=color ./ 255)
end

# PLOTTING SITES IN 3D (but complex)
function plotSiteComplex(
            site    :: S,
            radius  :: Real,
            color   :: Vector{<:Integer}
            ;
            detail  :: Int64 = 7,
            kwargs...
        ) where {L,S<:AbstractSite{L,3}}

    # get the coordinates
    x = point(site)[1]
    y = point(site)[2]
    z = point(site)[3]
    # Make data of sphere
    n = detail
    u = range(0,stop=2*pi,length=n);
    v = range(0,stop=pi,length=n);
    px = cos.(u) * sin.(v)';
    py = sin.(u) * sin.(v)';
    pz = ones(n) * cos.(v)';
    px = (px.*radius) .+ x
    py = (py.*radius) .+ y
    pz = (pz.*radius) .+ z
    # Plot the surface of the sphere
    PyPlot.surf(px,py,pz, rstride=1, cstride=1, linewidth=0, antialiased=true, color=color ./ 255)
end

# PLOTTING SITES IN 2D
function plotSite(
            site    :: S,
            radius  :: Real,
            color   :: Vector{<:Integer}
            ;
            site_labels :: Bool = true,
            site_label_fontsize :: Real = 12,
            site_label_offset :: Vector{<:Real} = [0.1, 0.0],
            kwargs...
        ) where {L,S<:AbstractSite{L,2}}

    # get the coordinates
    x = point(site)[1]
    y = point(site)[2]
    # scatter the point
    PyPlot.scatter([x,],[y,], color=color./255)
    # maybe annotate the label as text
    if site_labels
        PyPlot.text(x+site_label_offset[1],y+site_label_offset[2], string(label(site)), fontsize=site_label_fontsize, color=color./255)
    end
end

# PLOTTING SITES IN 3D
function plotSite(
            site    :: S,
            radius  :: Real,
            color   :: Vector{<:Integer}
            ;
            site_labels :: Bool = true,
            site_label_fontsize :: Real = 12,
            site_label_offset :: Vector{<:Real} = [0.1, 0.0, 0.0],
            kwargs...
        ) where {L,S<:AbstractSite{L,3}}

    # get the coordinates
    x = point(site)[1]
    y = point(site)[2]
    z = point(site)[3]
    # scatter the point
    PyPlot.scatter3D([x,],[y,],[z,], color=color./255)
    # maybe annotate the label as text
    if site_labels
        PyPlot.text3D(x+site_label_offset[1],y+site_label_offset[2],z+site_label_offset[3], string(label(site)), fontsize=site_label_fontsize, color=color./255)
    end
end

# plotting based on scattering

# 2D
function plotSites(
            sites   :: Vector{S},
            radius  :: Real,
            color   :: Vector{<:Integer}
            ;
            site_labels :: Bool = true,
            site_label_fontsize :: Real = 12,
            site_label_offset :: Vector{<:Real} = [0.1, 0.0],
            kwargs...
        ) where {L,S<:AbstractSite{L,2}}

    # scatter the points
    PyPlot.scatter([point(s)[1] for s in sites],[point(s)[2] for s in sites], color=color./255, zorder=10)
    # maybe annotate the label as text
    if site_labels
        for site in sites
            PyPlot.text(point(site)[1]+site_label_offset[1],point(site)[2]+site_label_offset[2], string(label(site)), fontsize=site_label_fontsize, color=color./255)
        end
    end
end

# 3D
function plotSites(
            sites   :: Vector{S},
            radius  :: Real,
            color   :: Vector{<:Integer}
            ;
            site_labels :: Bool = true,
            site_label_fontsize :: Real = 12,
            site_label_offset :: Vector{<:Real} = [0.1, 0.0, 0.0],
            kwargs...
        ) where {L,S<:AbstractSite{L,3}}

    # scatter the points
    PyPlot.scatter3D([point(s)[1] for s in sites],[point(s)[2] for s in sites],[point(s)[3] for s in sites], color=color./255)
    # maybe annotate the label as text
    if site_labels
        for site in sites
            PyPlot.text3D(point(site)[1]+site_label_offset[1],point(site)[2]+site_label_offset[2],point(site)[3]+site_label_offset[3], string(label(site)), fontsize=site_label_fontsize, color=color./255)
        end
    end
end

# 3D

# export the functions
export plotSite, plotSites, plotSiteComplex

################################################################################
#
#   PLOTTING BONDS
#
################################################################################

# PLOTTING BONDS IN 2D
function plotBond(
            site_from :: S,
            site_to   :: S,
            radius    :: Real,
            color     :: Vector{<:Integer}
            ;
            kwargs...
        ) where {L,S<:AbstractSite{L,2}}

    # get the coordinates
    x1 = point(site_from)[1]
    y1 = point(site_from)[2]
    # get the coordinates
    x2 = point(site_to)[1]
    y2 = point(site_to)[2]
    # plot a bond
    PyPlot.plot([x1,x2], [y1,y2], color=color./255)
end

# PLOTTING BONDS IN 3D
function plotBond(
            site_from :: S,
            site_to   :: S,
            radius    :: Real,
            color     :: Vector{<:Integer}
            ;
            kwargs...
        ) where {L,S<:AbstractSite{L,3}}

    # get the coordinates
    x1 = point(site_from)[1]
    y1 = point(site_from)[2]
    z1 = point(site_from)[3]
    # get the coordinates
    x2 = point(site_to)[1]
    y2 = point(site_to)[2]
    z2 = point(site_to)[3]
    # plot a bond
    PyPlot.plot3D([x1,x2], [y1,y2], [z1,z2], color=color./255)
end

# Plot bonds based on line collections
# https://gist.github.com/gizmaa/7214002#linecollection

# 2D
function plotBonds(
            bonds  :: Vector{B},
            sites  :: Vector{S},
            radius :: Real,
            color  :: Vector{<:Integer}
            ;
            kwargs...
        ) where {LS,S<:AbstractSite{LS,2},N,LB,B<:AbstractBond{LB,N}}

    # make a list of all lines
    lines = Vector{Tuple{Float64,Float64}}[
        Tuple{Float64,Float64}[(point(sites[from(b)])[1],point(sites[from(b)])[2]), (point(sites[to(b)])[1],point(sites[to(b)])[2])]
        for b in bonds
    ]
    # make the collection
    line_collection = matplotlib.collections.LineCollection(lines,colors=[color./255 for i in 1:length(lines)], label="$(label(bonds[1]))")
    # plot a bond
    PyPlot.gca().add_collection(line_collection)
end

# 3D NOT IMPLEMENTED YET (so tempoerary fix)
function plotBonds(
            bonds  :: Vector{B},
            sites  :: Vector{S},
            radius :: Real,
            color  :: Vector{<:Integer}
            ;
            kwargs...
        ) where {LS,S<:AbstractSite{LS,3},N,LB,B<:AbstractBond{LB,N}}

    # plot all bonds individually
    for b in bonds
        plotBond(
            sites[from(b)],
            sites[to(b)],
            radius,
            color;
            kwargs...
        )
    end
end

# export the functions
export plotBond, plotBonds


text:
include("utils.jl")
using MAT
using PyPlot
using ADCMEKit

## Load inital state from outside simulation Scycle
meta = matopen("data-scycle.mat")
data = read(meta, "d")
time = data["time"]
tau = data["tauP"]
v = data["slipVel"]
psi = data["psi"]
bd_right = data["momBal"]["bcR"]
bd_left = data["momBal"]["bcL"]
nt = 601
nx = 51
ny = 51
# eq_ind = collect(1:nt)[(v[1,:] .> 1e-1)[:]]
eq_ind = collect(190:200)

## model setup
m = 50
n = 50
h = 0.1

right = bcnode("right", m, n, h)
left = bcnode("left", m, n, h)
ind = [right;left]
ρ = 3
μ = 30
VPL = 1e-9 # plat loading velocity
f0 = 0.6

# Initial data
V0 = v[:, eq_ind[1]]
Ψ0 = psi[:, eq_ind[1]]
# Ψ0 = ones(n+1) * 1.01 * f0

η = constant(4.7434)
a = constant(ones(n+1) * 0.01)
b = constant(ones(n+1) * 0.02)
Dc = constant(0.03)
v0 = constant(1e-6)
f0 = constant(f0)
σn = constant(ones(n+1) * 50)
bd_left0 = constant(bd_left[:,eq_ind[1]])

function RangeVariable(a, b, n=nothing)
    if isnothing(n)
        (b-a)*sigmoid(Variable(0.0)) + a
    else
        (b-a)*sigmoid(Variable(zeros(n))) + a
    end
end

# η = RangeVariable(0.0,2.0) * 10.0
# a = RangeVariable(0.0, 2.0, n+1) * 0.02
# b = RangeVariable(0.0, 2.0, n+1) * 0.03
# Dc = RangeVariable(0.0,2.0)
# v0 = exp(-RangeVariable(0.0,8.0)*log(10))
# f0 = RangeVariable(0.0,2.0) * 0.8
# σn = RangeVariable(0.0, 2.0, n+1) * 30
# bd_left0 = RangeVariable(0.0, 2.0, n+1) * 0.3

# pl = placeholder(ones(1))
Δt = 0.01
t = collect(0:Δt:time[eq_ind[end]]-time[eq_ind[1]])
NT = length(t)
Δt = constant(ones(NT) .* Δt)

# @show NT
# NT = 5

# copy matrix to C++
ccall((:factorize_matrix, "./accel/build/libSaver"), Cvoid, (Cint, Cint, Cdouble, Ref{Cint}, Cint),
            Int32(m), Int32(n), h, Int32.(ind .- 1) , Int32(length(ind)))

function fast_solve(rhs)
    fast_solve_ = load_op_and_grad("./accel/build/libFastSolve","fast_solve")
    rhs = convert_to_tensor(rhs, dtype = Float64)
    fast_solve_(rhs)
end

K_ = compute_fem_stiffness_matrix1(diagm(0=>ones(2)), m, n, h)
K_[ind, :] .= 0.0
K_[ind, ind] = spdiagm(0=>ones(length(ind)))
K_ = dropzeros(K_)

## DEBUG
# u_bd = zeros(NT+1, n+1)
# for i = 1:NT+1
#     u_bd[i, :] = (i-1)*fix_Δt*VPL .+ bd_right[:, eq_ind[1]]
#     # u_bd[i, :] = bd_right[:, eq_ind[i]]
# end
# u_bd = constant(u_bd)

σzx_id = Int64[]
for i = 1:n
    elem = (i-1)*m + 1
    push!(σzx_id, (elem-1)*4 + 1)
end
push!(σzx_id, σzx_id[end]+2)

sess = Session();

# u0 = vector([left; right], [bd_left[:,eq_ind[1]]; bd_right[:,eq_ind[1]]], (m+1)*(n+1))
# U0 = K_\run(sess, u0)

# bd_left0 = Variable(bd_left[:,eq_ind[1]])
u0 = vector([left; right], [bd_left0; constant(bd_right[:,eq_ind[1]])], (m+1)*(n+1))
U0 = constant(K_)\u0

ψ_ref = constant(Array(psi[:,eq_ind]'))
bd_left_ref = constant(Array(bd_left[:,eq_ind]'))
tau_ref = constant(Array(tau[:,eq_ind]'))
v_ref = constant(Array(v[:,eq_ind]'))

# du is defined on the left side
function StateRate(Ψ, v)
    # v = abs(v)
    # s = v/(2v0)
    # o = exp(-Ψ/a)
    # f = a * (Ψ/a + log(s+sqrt(s^2+o^2)))
    f = a* tf.asinh(v/(2v0)*exp(Ψ/a))
    Δψ = -v/Dc*(f-f0+(b-a)*log(v/v0))
end

# https://acadpubl.eu/jsi/2015-101-5-6-7-8/2015-101-8/21/21.pdf
function rk3(i, u, psi, v_bd, Δt)
    ## --- RK Stage 1 ---
    du_bd_S1 = v_bd/2  # t

    # v_bd = v_ref[i] ## debug

    dpsi_S1 = StateRate(psi, v_bd) # t
    u_left_bd_S1 = u[left] + 0.5*Δt*du_bd_S1 # t+0.5Δt
    psi_S1 = psi + 0.5*Δt * dpsi_S1 # t+0.5Δt

    # psi_S1 = ψ_ref[i] ## debug
    # u_left_bd_S1 = bd_left_ref[i] ## debug

    u_right_bd_S1 = u[right] + 0.5*Δt*VPL # t + 0.5Δt
    u0 = vector([left;right], [u_left_bd_S1; u_right_bd_S1], (m+1)*(n+1))
    # u_S1 = K\u0 # t + 0.5Δt
    u_S1 = fast_solve(u0)

    ε = eval_strain_on_gauss_pts1(u_S1, m, n, h)
    σ = μ * ε
    τ = σ[σzx_id ,1]

    # τ = tau_ref[i] ## debug

    v_bd_S1 = compute_vel(a, v0, psi_S1, σn, τ, η)

    ## --- RK Stage 2 ---
    du_bd_S2 = v_bd_S1/2  # t + 0.5Δt

    # v_bd_S1 = v_ref[i] ## debug
    dpsi_S2 = StateRate(psi_S1, v_bd_S1) # t + 0.5Δt
    u_left_bd_S2 = u[left] + Δt * (-du_bd_S1 + 2*du_bd_S2) # t + Δt
    psi_S2 = psi + Δt * (-dpsi_S1 + 2*dpsi_S2) # t + Δt

    # psi_S2 = ψ_ref[i] ## debug
    # u_left_bd_S2 = bd_left_ref[i] ## debug

    u_right_bd_S2 = u[right] + Δt*VPL # t + Δt
    u0 = vector([left;right], [u_left_bd_S2; u_right_bd_S2], (m+1)*(n+1))# t + Δt
    # u_S2 = K\u0 # t + Δt
    u_S2 = fast_solve(u0)

    ε = eval_strain_on_gauss_pts1(u_S2, m, n, h)
    σ = μ * ε
    τ = σ[σzx_id,1]

    ## debug
    # τ = tau_ref[i]
    # op = tf.print(τ)
    # i = bind(i, op)

    v_bd_S2 = compute_vel(a, v0, psi_S2, σn, τ, η) # t+ Δt + δΔt

    ## --- RK Stage 3 ---
    du_bd_S3 = v_bd_S2/2 # t + Δt

    # v_bd_S2 = v_ref[i] ## debug
    dpsi_S3  = StateRate(psi_S2, v_bd_S2) # t + Δt

    u_left_bd_S3 = u[left] + Δt/6*(du_bd_S1 + 4*du_bd_S2 + du_bd_S3); ## t + Δt
    psi_S3 = psi + Δt/6*(dpsi_S1 + 4*dpsi_S2 + dpsi_S3); ## t + Δt

    u_right_bd_S3 = u[right] + Δt*VPL # t + Δt
    u0 = vector([left;right], [u_left_bd_S3; u_right_bd_S3], (m+1)*(n+1))# t + Δt
    # u_S3 = K\u0 # t + Δt
    u_S3 = fast_solve(u0)

    ε = eval_strain_on_gauss_pts1(u_S3, m, n, h)
    σ = μ * ε
    τ = σ[σzx_id,1]

    # τ = tau_ref[i+1] ## debug
    # psi_S3 = ψ_ref[i+1] ## debug

    v_bd_S3 = compute_vel(a, v0, psi_S3, σn, τ, η) # t+ Δt + δΔt

    # psi_S3 = ψ_ref[i+1] ## debug
    # u_left_bd = bd_left_ref[i+1]  ## debug

    return u_S3, psi_S3, v_bd_S3
end

function simulate()

    function body(i, ta_u, ta_Ψ, ta_v_bd)
        u = read(ta_u, i)
        Ψ = read(ta_Ψ, i)
        v_bd = read(ta_v_bd, i)
        u, Ψ, v_bd = rk3(i, u, Ψ, v_bd, Δt[i])
        ta_u = write(ta_u, i+1, u)
        ta_Ψ = write(ta_Ψ, i+1, Ψ)
        # ta_Ψ = write(ta_Ψ, i+1, ψ_ref[i+1]) ## debug
        ta_v_bd = write(ta_v_bd, i+1, v_bd)
        return i+1, ta_u, ta_Ψ, ta_v_bd
    end

    function condition(i, ta_u, ta_Ψ, ta_v_bd)
        i<=NT
    end

    ta_u, ta_Ψ, ta_v_bd = TensorArray(NT+1), TensorArray(NT+1), TensorArray(NT+1)
    i = constant(1, dtype=Int32)
    ta_u = write(ta_u, 1, convert_to_tensor(U0))
    # ta_u = write(ta_u, 1, constant(U0))
    ta_Ψ = write(ta_Ψ, 1, constant(Ψ0))
    ta_v_bd = write(ta_v_bd, 1, constant(V0))

    _, u, Ψ, v_bd = while_loop(condition,body, [i, ta_u, ta_Ψ, ta_v_bd])
    u, Ψ, v_bd = stack(u) , stack(Ψ),  stack(v_bd)
    u = set_shape(u, (NT+1, size(u,2)))
    return u, Ψ, v_bd
end

## debug
# using PyCall
# py"""
# import traceback
# try:
#     # Your codes here
#     $sess.run($u)
# except Exception:
#     print(traceback.format_exc())
# """

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages