text:
using DiffEqJump, DiffEqBase, OrdinaryDiffEq
using Test, Graphs, LinearAlgebra

reltol       = 0.05
Nsims = 10^4

dim = 1
linear_size = 5
dims = Tuple(repeat([linear_size], dim))
num_nodes = prod(dims)
center_site = trunc(Int,(linear_size^dim + 1)/2)
u0 = zeros(Int, 1, num_nodes)
end_time = 100.0
diffusivity = 1.0
death_rate = 0.1

num_species = 1
reactstoch = [Pair{Int64, Int64}[], [1 => 1]]
netstoch = [[1 => 1],[1 => -1]]
uniform_rates = ones(2, num_nodes)
uniform_rates[2,:] *= death_rate
non_uniform_rates = zeros(2, num_nodes)
non_uniform_rates[:,center_site] = uniform_rates[:,center_site]

# DiscreteProblem setup
tspan = (0.0, end_time)
prob = DiscreteProblem(u0, tspan, uniform_rates)

# spatial system
grid = Graphs.grid(dims)
hopping_constants = [diffusivity for i in u0]

# majumps
uniform_majumps_1 = SpatialMassActionJump(uniform_rates[:,1], reactstoch, netstoch)
uniform_majumps_2 = SpatialMassActionJump(uniform_rates, reactstoch, netstoch)
uniform_majumps_3 = SpatialMassActionJump([1.], reshape(uniform_rates[2,:], 1, num_nodes), reactstoch, netstoch) # hybrid
uniform_majumps_4 = SpatialMassActionJump(MassActionJump(uniform_rates[:,1], reactstoch, netstoch))
uniform_majumps = [uniform_majumps_1, uniform_majumps_2, uniform_majumps_3, uniform_majumps_4]

non_uniform_majumps_1 = SpatialMassActionJump(non_uniform_rates, reactstoch, netstoch) # reactions are zero outside of center site
non_uniform_majumps_2 = SpatialMassActionJump([1.], reshape(non_uniform_rates[2,:], 1, num_nodes), reactstoch,netstoch) # birth everywhere, death only at center site
non_uniform_majumps_3 = SpatialMassActionJump([1. 0. 0. 0. 0.; 0. 0. 0. 0. death_rate], reactstoch,netstoch) # birth on the left, death on the right
non_uniform_majumps = [non_uniform_majumps_1, non_uniform_majumps_2, non_uniform_majumps_3]

# put together the JumpProblem's
uniform_jump_problems = JumpProblem[JumpProblem(prob, NSM(), majump, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false)) for majump in uniform_majumps]
# flattenned
append!(uniform_jump_problems, JumpProblem[JumpProblem(prob, NRM(), majump, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false)) for majump in uniform_majumps])

# non-uniform
non_uniform_jump_problems = JumpProblem[JumpProblem(prob, NSM(), majump, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false)) for majump in non_uniform_majumps]

# testing
function get_mean_end_state(jump_prob, Nsims)
    end_state = zeros(size(jump_prob.prob.u0))
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        end_state .+= sol.u[end]
    end
    end_state/Nsims
end

function discrete_laplacian_from_spatial_system(spatial_system, hopping_rate)
    sites = 1:num_sites(spatial_system)
    laplacian = zeros(length(sites), length(sites))
    for site in sites
        laplacian[site,site] = -outdegree(spatial_system, site)
        for nb in neighbors(spatial_system, site)
            laplacian[site, nb] = 1
        end
    end
    laplacian .*= hopping_rate
    laplacian
end
L = discrete_laplacian_from_spatial_system(grid, diffusivity)

# birth and death everywhere
f(u,p,t) = L*u - death_rate*u + uniform_rates[1,:]
ode_prob = ODEProblem(f, zeros(num_nodes), tspan)
sol = solve(ode_prob, Tsit5())

for spatial_jump_prob in uniform_jump_problems
    solution = solve(spatial_jump_prob, SSAStepper())
    mean_end_state = get_mean_end_state(spatial_jump_prob, Nsims)
    mean_end_state = reshape(mean_end_state, num_nodes)
    diff = mean_end_state - sol.u[end]
    for (i,d) in enumerate(diff)
        @test abs(d) < reltol*sol.u[end][i]
    end
end

# birth and death zero outside of center site
f(u,p,t) = L*u - diagm([0., 0., death_rate, 0., 0.])*u + [0., 0., 1., 0., 0.]
ode_prob = ODEProblem(f, zeros(num_nodes), tspan)
sol = solve(ode_prob, Tsit5())

solution = solve(non_uniform_jump_problems[1], SSAStepper())
mean_end_state = get_mean_end_state(non_uniform_jump_problems[1], Nsims)
mean_end_state = reshape(mean_end_state, num_nodes)
diff = mean_end_state - sol.u[end]
for (i,d) in enumerate(diff)
    @test abs(d) < reltol*sol.u[end][i]
end

# birth everywhere, death only at center site
f(u,p,t) = L*u - diagm([0., 0., death_rate, 0., 0.])*u + ones(num_nodes)
ode_prob = ODEProblem(f, zeros(num_nodes), tspan)
sol = solve(ode_prob, Tsit5())

solution = solve(non_uniform_jump_problems[2], SSAStepper())
mean_end_state = get_mean_end_state(non_uniform_jump_problems[2], Nsims)
mean_end_state = reshape(mean_end_state, num_nodes)
diff = mean_end_state - sol.u[end]
for (i,d) in enumerate(diff)
    @test abs(d) < reltol*sol.u[end][i]
end

# birth on left end, death on right end
f(u,p,t) = L*u - diagm([0., 0., 0., 0., death_rate])*u + [1., 0., 0., 0., 0.]
ode_prob = ODEProblem(f, zeros(num_nodes), tspan)
sol = solve(ode_prob, Tsit5())

solution = solve(non_uniform_jump_problems[3], SSAStepper())
mean_end_state = get_mean_end_state(non_uniform_jump_problems[3], Nsims)
mean_end_state = reshape(mean_end_state, num_nodes)
diff = mean_end_state - sol.u[end]
for (i,d) in enumerate(diff)
    @test abs(d) < reltol*sol.u[end][i]
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: pinpoint_sms_voice
using AWS.Compat
using AWS.UUIDs

"""
    create_configuration_set()
    create_configuration_set(params::Dict{String,<:Any})

Create a new configuration set. After you create the configuration set, you can add one or
more event destinations to it.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ConfigurationSetName"`: The name that you want to give the configuration set.
"""
function create_configuration_set(; aws_config::AbstractAWSConfig=global_aws_config())
    return pinpoint_sms_voice(
        "POST",
        "/v1/sms-voice/configuration-sets";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_configuration_set(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return pinpoint_sms_voice(
        "POST",
        "/v1/sms-voice/configuration-sets",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_configuration_set_event_destination(configuration_set_name)
    create_configuration_set_event_destination(configuration_set_name, params::Dict{String,<:Any})

Create a new event destination in a configuration set.

# Arguments
- `configuration_set_name`: ConfigurationSetName

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"EventDestination"`:
- `"EventDestinationName"`: A name that identifies the event destination.
"""
function create_configuration_set_event_destination(
    ConfigurationSetName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return pinpoint_sms_voice(
        "POST",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_configuration_set_event_destination(
    ConfigurationSetName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return pinpoint_sms_voice(
        "POST",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_configuration_set(configuration_set_name)
    delete_configuration_set(configuration_set_name, params::Dict{String,<:Any})

Deletes an existing configuration set.

# Arguments
- `configuration_set_name`: ConfigurationSetName

"""
function delete_configuration_set(
    ConfigurationSetName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return pinpoint_sms_voice(
        "DELETE",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_configuration_set(
    ConfigurationSetName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return pinpoint_sms_voice(
        "DELETE",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_configuration_set_event_destination(configuration_set_name, event_destination_name)
    delete_configuration_set_event_destination(configuration_set_name, event_destination_name, params::Dict{String,<:Any})

Deletes an event destination in a configuration set.

# Arguments
- `configuration_set_name`: ConfigurationSetName
- `event_destination_name`: EventDestinationName

"""
function delete_configuration_set_event_destination(
    ConfigurationSetName,
    EventDestinationName;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return pinpoint_sms_voice(
        "DELETE",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations/$(EventDestinationName)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_configuration_set_event_destination(
    ConfigurationSetName,
    EventDestinationName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return pinpoint_sms_voice(
        "DELETE",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations/$(EventDestinationName)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_configuration_set_event_destinations(configuration_set_name)
    get_configuration_set_event_destinations(configuration_set_name, params::Dict{String,<:Any})

Obtain information about an event destination, including the types of events it reports,
the Amazon Resource Name (ARN) of the destination, and the name of the event destination.

# Arguments
- `configuration_set_name`: ConfigurationSetName

"""
function get_configuration_set_event_destinations(
    ConfigurationSetName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return pinpoint_sms_voice(
        "GET",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_configuration_set_event_destinations(
    ConfigurationSetName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return pinpoint_sms_voice(
        "GET",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_configuration_sets()
    list_configuration_sets(params::Dict{String,<:Any})

List all of the configuration sets associated with your Amazon Pinpoint account in the
current region.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"NextToken"`: A token returned from a previous call to the API that indicates the
  position in the list of results.
- `"PageSize"`: Used to specify the number of items that should be returned in the response.
"""
function list_configuration_sets(; aws_config::AbstractAWSConfig=global_aws_config())
    return pinpoint_sms_voice(
        "GET",
        "/v1/sms-voice/configuration-sets";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_configuration_sets(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return pinpoint_sms_voice(
        "GET",
        "/v1/sms-voice/configuration-sets",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    send_voice_message()
    send_voice_message(params::Dict{String,<:Any})

Create a new voice message and send it to a recipient's phone number.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CallerId"`: The phone number that appears on recipients' devices when they receive the
  message.
- `"ConfigurationSetName"`: The name of the configuration set that you want to use to send
  the message.
- `"Content"`:
- `"DestinationPhoneNumber"`: The phone number that you want to send the voice message to.
- `"OriginationPhoneNumber"`: The phone number that Amazon Pinpoint should use to send the
  voice message. This isn't necessarily the phone number that appears on recipients' devices
  when they receive the message, because you can specify a CallerId parameter in the request.
"""
function send_voice_message(; aws_config::AbstractAWSConfig=global_aws_config())
    return pinpoint_sms_voice(
        "POST",
        "/v1/sms-voice/voice/message";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function send_voice_message(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return pinpoint_sms_voice(
        "POST",
        "/v1/sms-voice/voice/message",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_configuration_set_event_destination(configuration_set_name, event_destination_name)
    update_configuration_set_event_destination(configuration_set_name, event_destination_name, params::Dict{String,<:Any})

Update an event destination in a configuration set. An event destination is a location that
you publish information about your voice calls to. For example, you can log an event to an
Amazon CloudWatch destination when a call fails.

# Arguments
- `configuration_set_name`: ConfigurationSetName
- `event_destination_name`: EventDestinationName

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"EventDestination"`:
"""
function update_configuration_set_event_destination(
    ConfigurationSetName,
    EventDestinationName;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return pinpoint_sms_voice(
        "PUT",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations/$(EventDestinationName)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_configuration_set_event_destination(
    ConfigurationSetName,
    EventDestinationName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return pinpoint_sms_voice(
        "PUT",
        "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations/$(EventDestinationName)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

ext: jl


text:
using Base.Test
using PropertyGraph
using UUID

# Uncomment this include when running this test directly....
# this include is not required when running vi runtests.jl
# include("testdata.jl")

# run the test within a function scope to avoid interference with other tests
function compositequerytest()
	# build a test graph that will be used to test graph queries
	g = buildstatetestgraph()

	# start with Queensland and Victoria
	q = vertices(g, v-> get(v,"name") == "Queensland" || get(v,"name") == "Victoria")
	# follow outgoung edges of type AdjacentTo
	q = outgoing(q, e->e.typelabel == "IsAdjacent")
	# and include the vertices the edges lead to
	q = head(q)
	# reduce the set to the distinct set
	q = distinct(q)

	# count adjacent states
	countadjacentstates = count(q)
	# expect 3 adjacent states (new south wales, northern territory, south australia)
	@test countadjacentstates == 3

	# and take the sum of the population
	sumadjacentpopulations = sum(q, v->get(v, "population"))
	# expect the sum of the populations from (new south wales, northern territory, south australia)
	@test sumadjacentpopulations == 9355700

	# retrieve the minimum of the population
	minimumadjacentpopulation = minimum(q, v->get(v, "population"))
	# expect the minimum of the adjacent populations to be the population of northern territory
	@test minimumadjacentpopulation == 241800

	# retrieve the maximum of the population
	maximumadjacentpopulation = maximum(q, v->get(v, "population"))
	# expect the maximum of the adjacent populations to be the population of new south wales
	@test maximumadjacentpopulation == 7439200

	# calculate the average of the adjacent populations
	averageadjacentpopulation = average(q, v->get(v, "population"))
	# expect the maximum of the adjacent populations to be the population of new south wales
	@test ifloor(averageadjacentpopulation) == 3118566

	# test short form queries
	countadjacentstates = query(g,
			  (vertices, v-> get(v,"name") == "Queensland" || get(v,"name") == "Victoria"),
			  (outgoing, e->e.typelabel == "IsAdjacent"),
			  head,
			  distinct,
			  count
			 )
	@test countadjacentstates == 3

	averageadjacentpopulation = query(g,
				 (vertices, v-> get(v,"name") == "Queensland1" || get(v,"name") == "Victoria1"),
				 (outgoing, e->e.typelabel == "IsAdjacent"),
				 head,
				 distinct,
				 (average, v->get(v, "population"))
			)
	@test averageadjacentpopulation == UnspecifiedValue

	#test building the query in parts
	q = query(g,
				(vertices, v-> get(v,"name") == "Queensland" || get(v,"name") == "Victoria"),
			 	(outgoing, e->e.typelabel == "IsAdjacent"),
			 )
	# build the rest of the query
	averageadjacentpopulation = query(q,
				 head,
				 distinct,
				 (average, v->get(v, "population"))
			 )
	@test ifloor(averageadjacentpopulation) == 3118566

	# get the state capitals whose population is > 10% of their nations population
	# this query has the following structure
	#       start with a graph
	#       select vertices of type "Country"
	#       store the countries population (so that it will be carried forward with graph traversals for later comparison)
	#       follow the incoming edges of type IsIn to move to States of each Country
	#       the tail of each edge is a State
	#       follow the incoming edges of type IsCaptitalOf
	#       the tail of these edges are capital cities, filter to only those with a population > 10% of stored country population
	q = query(g,
				(vertices, v-> v.typelabel == "Country"),
			 	(store, v->@compat Dict{String, Any}("countrypopulation"=>get(v,"population",0),"countryname"=>get(v,"name"))),
			  	(incoming, e->e.typelabel == "IsIn"),
				tail,
			  	(incoming, e->e.typelabel == "IsCapitalOf"),
				(tail, v->get(v,"population",0) > getstored(v,"countrypopulation",0) * 0.10)
			 )
	countcities = query(q, count)
	@test countcities == 2 # in this test, 2 cities are expected

	# select a flat set of results in the form:  cityname, citypopulation, countryname, countrypopulation
	selectresults = query(q, (select, s->(get(s,"name"),
								 get(s,"population"),
								 getstored(s,"countryname"),
								 getstored(s,"countrypopulation"))))

	@test length(selectresults) == countcities

	# test the data expected was returned
	sydneymatched = false
	melbournematched = false

	for r in selectresults
		if r[1] == "Sydney"
			sydneymatched = true
			@test r[1] == "Sydney"
			@test r[2] == 4667283
			@test r[3] == "Australia"
			@test r[4] == 23235800
		end
		if r[1] == "Melbourne"
			melbournematched = true
			@test r[1] == "Melbourne"
			@test r[2] == 4246345
			@test r[3] == "Australia"
			@test r[4] == 23235800
		end
	end

	@test sydneymatched
	@test melbournematched

	# find the States whose capitals have more than half the State's population
	# this can be achieved in several ways... however for this stest Store and Move is being targetted for testing
	# once the cities are identified, movetovertex is used to move back to previously stored State vertices
	populationquery = query(g,
				(vertices, v-> v.typelabel == "State"),
			 	(store, v->@compat Dict{String, Any}("state"=>v.item, "statepopulation"=>get(v, "population"))),
			  	(incoming, e->e.typelabel == "IsCapitalOf"),
				(tail, v->get(v,"population",0) > (getstored(v,"statepopulation",0) * 0.5)),
			  	(movetovertices, v->getstored(v,"state")) # move back to State
			 )

	sumpopulation = query(populationquery, (sum, v->get(v,"population")))
	@test sumpopulation == 17931600

	selectpopulation = query(populationquery, (select, s->(get(s,"name"), get(s,"population"))))
	for r in selectpopulation
		println("State: ", r[1], " Population: ", r[2])
	end

	# repeat the tests above, but this time move to edges instead of vertices
	sumpopulation = query(g,
				(vertices, v-> v.typelabel == "State"),
			 	(store, v->@compat Dict{String, Any}("statepopulation"=>get(v, "population"))),
			  	(incoming, e->e.typelabel == "IsCapitalOf"),
				(store, e->@compat Dict{String, Any}("capitaledge"=>e.item)),
			  	(tail, v->get(v,"population",0) > getstored(v,"statepopulation",0) * 0.5),
			  	(movetoedges, v->getstored(v,"capitaledge")), # move back to State via the stored edge
			  	head,
				(sum, v->get(v,"population"))
			 )
	@test sumpopulation == 17931600

	# sum populations, grouped by state or territory type
	groupedsum = query(g,
					   (vertices, v-> v.typelabel == "State" || v.typelabel == "Territory"),
					   (group, v-> v.typelabel),
					   (sum, v->get(v,"population")),
					   )
	# test grouped sum
	@test groupedsum["Territory"] == 624700

	# sum populations of capitals of states or territories, grouped by state or territory type
	groupedsum = query(g,
					   (vertices, v-> v.typelabel == "State" || v.typelabel == "Territory"),
					   (group, v-> v.typelabel),
					   (incoming, e->e.typelabel == "IsCapitalOf"),
					   tail,
					   (sum, v->get(v,"population"))
					   )
	@test groupedsum["Territory"] == 514578

	# multi-level group
	# sum the population of adjacent states for each state by adjacency direction
	groupedsum = query(g,
					   (vertices, v-> v.typelabel == "State" || v.typelabel == "Territory"),
					   (group, v-> get(v,"name")),
					   (outgoing, e->e.typelabel == "IsAdjacent"),
					   (group, e-> get(e,"direction","none")),
					   head,
					   (sum, v->get(v,"population"))
					   )

	# test the sum of populations of States directly West of Queensland
	@test groupedsum["Queensland"]["West"] == 1916500
end

compositequerytest()


text:
module Data_module
    const CHEMICAL_NUMBER = "chemical.number"

    const EPS_DEFAULT = 1.0E-15

    const GAUSS_LEGENDRE_INTEG_DEFAULT = 2

    const GRID_NUM_DEFAULT = 20000

    const ITERATION_CRITERION_DEFAULT = 1.0E-13

    const ITERATION_MAXITER_DEFAULT = 1000

    const ITERATION_MIXING_WEIGHT_DEFAULT = 0.1

    const XMAX_DEFAULT = 100.0

    const XMIN_DEFAULT = 1.0E-5

    const X_MATCHING_POINT_DEFAULT = 11.0

    mutable struct Data_val
        usethread::Bool
        eps::Union{Float64, Nothing}
        gauss_legendre_integ_num::Union{Int32, Nothing}
        grid_num::Union{Int32, Nothing}
        iteration_criterion::Union{Float64, Nothing}
        iteration_maxiter::Union{Int32, Nothing}
        iteration_mixing_weight::Union{Float64, Nothing}
        xmax::Union{Float64, Nothing}
        xmin::Union{Float64, Nothing}
        x_matching_point::Union{Float64, Nothing}
        Z::Union{Float64, Nothing}
    end
end


text:
@testset "ParameterTable - RAMMatrices conversion" begin
    partable = ParameterTable(ram_matrices)
    @test ram_matrices == RAMMatrices(partable)
end

@test get_identifier_indices([:x2, :x10, :x28], model_ml) == [2, 10, 28]

@testset "get_identifier_indices" begin
    pars = [:θ_1, :θ_7, :θ_21]
    @test get_identifier_indices(pars, model_ml) == get_identifier_indices(pars, partable)
    @test get_identifier_indices(pars, model_ml) == get_identifier_indices(pars, RAMMatrices(partable))
end

graph = @StenoGraph begin
    # measurement model
    visual  → fixed(1.0, 1.0)*x1 + fixed(0.5,     0.5)*x2 + fixed(0.6, 0.8)*x3
    textual → fixed(1.0, 1.0)*x4 +                     x5 + label(:a₁, :a₂)*x6
    speed   → fixed(1.0, 1.0)*x7 + fixed(1.0,     NaN)*x8 + label(:λ₉, :λ₉)*x9
    # variances and covariances
    _(observed_vars) ↔ _(observed_vars)
    _(latent_vars)   ↔ _(latent_vars)
    visual ↔ textual + speed
    textual ↔ speed
end


text:
using InvariantMeasures
using ValidatedNumerics

using Plots
using LaTeXStrings

D = Mod1Dynamic(x -> 4*x + 0.01*InvariantMeasures.sinpi(8*x))
#D = Mod1Dynamic(x -> 16*x + 0.01*InvariantMeasures.sinpi(32*x))
#D = Mod1Dynamic(x->2*x+0.5*x*(1-x))
# D = PwMap(
# [x -> 17x/5,
#  x -> 34(x-5//17)^2/25 + 3(x-5//17),
#  x -> 34(x-10//17)^2/25 + 3(x-10//17),
#  x -> 17(x-15//17)/5
# ],
# [0, @interval(5/17), @interval(10/17), @interval(15/17), 1]
# )

function plot_norms(D; args...)
    dims = 2 .^ (1:12)
    num_norms = 30

    norms = fill(NaN, (num_norms, length(dims)))

    for (i, d) in enumerate(dims)
        B = Ulam(d)
        Q = DiscretizedOperator(B, D)
        norms[:,i] = norms_of_powers(weak_norm(B), num_norms, Q, integral_covector(B))
    end

    pgfplotsx()
    plot(norms,
        label = L"n = " .* string.(dims'),
        yscale= :log10,
        legend = :topright,
    #    ylims = (1e-6, 1)
        xlabel = L"k",
        ylabel = L"computational bounds to $\|Q^k|_U\|$";
        args...
    )
end

p1 = plot_norms(Mod1Dynamic(x -> 4*x + 0.01*InvariantMeasures.sinpi(8*x)),
                title = L"4x + 0.01\sin(8 \pi x)", legend=:none)

p2 = plot_norms(Mod1Dynamic(x->2*x+0.5*x*(1-x)),
                title = "Lanford map", legend = :bottomleft)

p = plot(p1, p2)

# savefig(p, "norms_vs_dimension_plot.tikz")


text:
using GalacticOptim, Optim, Test

rosenbrock(x, p) =  (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p  = [1.0, 100.0]

l1 = rosenbrock(x0, _p)
prob = OptimizationProblem(rosenbrock, x0, p=_p)
sol = solve(prob, SimulatedAnnealing())
@test 10*sol.minimum < l1

prob = OptimizationProblem(rosenbrock, x0, p=_p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, SAMIN())
@test 10*sol.minimum < l1

rosenbrock(x, p=nothing) =  (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2

l1 = rosenbrock(x0)
prob = OptimizationProblem(rosenbrock, x0)
sol = solve(prob, NelderMead())
@test 10*sol.minimum < l1

optprob = OptimizationFunction(rosenbrock, x0, GalacticOptim.AutoZygote())

prob = OptimizationProblem(optprob, x0)
sol = solve(prob, BFGS())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0)
sol = solve(prob, Newton())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0)
sol = solve(prob, Optim.KrylovTrustRegion())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, Fminbox())
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, SAMIN())
@test 10*sol.minimum < l1

using NLopt
prob = OptimizationProblem(optprob, x0)
sol = solve(prob, Opt(:LN_BOBYQA, 2))
@test 10*sol.minimum < l1

sol = solve(prob, Opt(:LD_LBFGS, 2))
@test 10*sol.minimum < l1

prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, Opt(:G_MLSL_LDS, 2), nstart=5, local_method = Opt(:LD_LBFGS, 2))
@test 10*sol.minimum < l1

# using MultistartOptimization
# sol = solve(prob, MultistartOptimization.TikTak(100), local_method = NLopt.LD_LBFGS)
# @test 10*sol.minimum < l1

# using QuadDIRECT
# sol = solve(prob, QuadDirect(); splits = ([-0.5, 0.0, 0.5],[-0.5, 0.0, 0.5]))
# @test 10*sol.minimum < l1

using Evolutionary
sol = solve(prob, CMAES(μ =40 , λ = 100),abstol=1e-15)
@test 10*sol.minimum < l1

using BlackBoxOptim
prob = GalacticOptim.OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
sol = solve(prob, BBO())
@test 10*sol.minimum < l1


text:
# Bui Quang Bao
# 19120454

# Data Mining - Lab02 - Frequent Itemset Mining

# How to use:
# 1. Change some settings below
# 2. Open command line, cd to directory has "apriori.jl"
# 3. Run: julia apriori.jl

# Settings ===================================
input_file_name = "data/test.txt"             # Input file
min_support_count = 3                         # Minimum support count for Apriori Algorithm
min_confidence = 0.8                          # Only choose rules have confidence >= min_confidence
                                              # To print all Association Rules: min_confidence = 0
# ============================================

# Just un-comment 1 group of codes below for quick settings:

# input_file_name = "data/mushrooms.txt"
# min_support_count = 4000
# min_confidence = 0.97

# input_file_name = "data/chess.txt"
# min_support_count = 3000
# min_confidence = 0.98

# input_file_name = "data/foodmartFIM.txt"      # Not recommended, take time to finish
# min_support_count = 15
# min_confidence = 0

# input_file_name = "data/retail.txt"           # Not recommended, take A LOT OF time to finish
# min_support_count = 10000
# min_confidence = 0.5

# ============================================

# Return True if set 1 is a subset of set 2, False if not
function isSubSet(l1, l2)
    for i in l1
        if !(i in l2)
            return false
        end
    end
    return true
end

# Return True if set 1 is equal to set 2, False if not
function isEqualSet(l1, l2)
    if length(l1) != length(l2)
        return false
    end
    for i in l1
        if !(i in l2)
            return false
        end
    end
    return true
end

# Sub-function for Apriori Algorithm
# Input: data list
# Return: L1
function apriori_data_to_L1(data, min_supp)
    items_list = []
    # Get all items and count
    for t in data
        for item in t
            existed = false
            for e in items_list
                if item == e[1][1]
                    existed = true
                    e[2] += 1
                end
            end
            if !existed
                push!(items_list, [[item], 1])
            end
        end
    end
    # Only accept items with count >= min_supp
    result_list = []
    for e in items_list
        if e[2] >= min_supp
            push!(result_list, e)
        end
    end
    return result_list
end
# Sub-function for Apriori Algorithm
# Input: L[k-1]
# Return: Items of C[k]
function apriori_L_to_C_items(L)
    print("L -> C's items")
    items_list = []
    result_list = []
    # Create an unique items list
    for e in L
        for i in e[1]
            if !(i in items_list)
                push!(items_list, i)
            end
        end
    end
    # Create combined list
    for e in L
        for i in items_list
            temp = []
            if !(i in e[1])
                for j in e[1]
                    push!(temp, j)
                end
                push!(temp, i)
                # temp.sort()
                sort!(temp)
                if !(temp in result_list)
                    push!(result_list, temp)
                end
            end
        end
    end
    return result_list
end
# Sub-function for Apriori Algorithm
# Input: Items of C[k]
# Return: C[k]
function apriori_C_items_to_C(data, C_items)
    print(" | C's items -> C")
    C = [[e, 0] for e in C_items]
    # Start counting
    for t in data
        for c in C
            if isSubSet(c[1],t)
                c[2] += 1
            end
        end
    end
    return C
end
# Sub-function for Apriori Algorithm
# Input: C[k]
# Return: L[k]
function apriori_C_to_L(C, min_supp)
    println(" | C -> L")
    L = []
    for c in C
        if c[2] >= min_supp
            push!(L, c)
        end
    end
    return L
end
# Sub-function for Apriori Algorithm
# Just print the apriori result to the console
function apriori_print(L)
    println("============================ Result ===============================")
    println("Itemset -> Support Count")
    for l in L
        for e in l
            ls_print(e[1])
            print(" -> ")
            println(e[2])
        end
    end
    println("===================================================================")
end

# Just to print array beautifully
function ls_print(ls)
    print("[")
    first_ele = true
    for e in ls
        if first_ele
            print("$e")
        else
            print(", $e")
        end
        first_ele = false
    end
    print("]")
end

# ======================== Apriori Algorithm =========================
function apriori(data, min_support_count)
    L = []
    # Initialize L with L1
    push!(L, apriori_data_to_L1(data, min_support_count))
    k = 2
    # Find all L[k] from L[k-1] until L[k-1] is empty
    while !isempty(L[k-1])
        print("> Processing: L[$(k-1)] to L[$(k)]: ")
        # Find L[k] from L[k-1]
        # Use all the sub-function for Apriori
        Lk = apriori_C_to_L(
            apriori_C_items_to_C(
                data, apriori_L_to_C_items(L[k-1])
            ), min_support_count
        )
        # Add L[k] to L
        push!(L, Lk)
        k += 1
    end
    # Just remove the last is an empty list
    pop!(L)
    # Just print the apriori result to the console
    apriori_print(L)
    return L
end
# ====================================================================

# Return the confidence of a rule: KB -> alpha
function apriori_confidence(apriori_result, KB, alpha)
    itemset = []
    for e in KB
        push!(itemset, e)
    end
    for e in alpha
        if !(e in itemset)
            push!(itemset, e)
        end
    end

    sup_count_itemset = 0
    sup_count_KB = -1

    for l in apriori_result
        for i in l
            if isEqualSet(itemset, i[1])
                sup_count_itemset = i[2]
            end
            if isEqualSet(KB, i[1])
                sup_count_KB = i[2]
            end
        end
    end

    # print("Confidence($KB => $alpha) = ")
    # print(sup_count_itemset/sup_count_KB)
    # println(" - $sup_count_itemset / $sup_count_KB")
    # println("===================================================================")

    return sup_count_itemset/sup_count_KB
end

# Return powerset of a set
# Example: [1, 2, 3]
# Return: [1], [2], [1, 2], [3], [1, 3], [2, 3], [1, 2, 3]
function powerset(x::Vector{T}) where T
    result = Vector{T}[[]]
    for e in x, j in eachindex(result)
        push!(result, [result[j] ; e])
    end
    return result
end

# Association Rules Generation
# Print all rules (from an itemset) with confidence >= min_confidence
function rules_gen(apriori_result, itemset, min_confidence)
    if length(itemset) != 0
        n = length(itemset)
        rules_right = []
        rules_left = []
        for powerset in powerset(itemset)
            if (length(powerset) != n) && (length(powerset) != 0)
                push!(rules_left, powerset)
            end
        end
        for e in rules_left
            push!(rules_right,symdiff(itemset, e))
        end
        # Print to console
        for i = 1:length(rules_left)
            conf = apriori_confidence(apriori_result, rules_left[i], rules_right[i])
            if conf >= min_confidence
                # print("Confidence($(rules_left[i]) => $(rules_right[i])) = ")
                print("Confidence(")
                ls_print(rules_left[i])
                print(" => ")
                ls_print(rules_right[i])
                print(") = ")
                println(round(conf; digits = 2))
            end
        end
    end
end

# Just apply rules_gen for all frequent itemset
function rules_gen_for_all_itemset(apriori_result, min_confidence)
    println("All Association Rules with confidence >= $min_confidence:")
    for i = 2:length(apriori_result)
        for e in apriori_result[i]
            rules_gen(apriori_result, e[1], min_confidence)
        end
    end
end

# Main
open(input_file_name) do f_input
    println("====================== Apriori Processing =========================")
    lines = []
    temp = []
    data = []
    # Read input file
    while ! eof(f_input)
       s = rstrip(readline(f_input))
       push!(lines, s)
    end
    println("> Read file successfully")
    # Turn data into data list
    for line in lines
        push!(temp, split(line, " "))
    end
    for t in temp
        push!(data, [parse(Int, x) for x in t])
    end
    println("> Turn data into proper format successfully")
    # ===========================================================================
    # Apriori
    apriori_result = apriori(data, min_support_count)
    # Rule Generation
    rules_gen_for_all_itemset(apriori_result, min_confidence)
    println("===================================================================")
end

max_stars_repo_path: CSC14004 - Data Mining and Application/Lab 02 - Frequent Itemset Mining - Apriori and FP-Tree (Julia)/Source/apriori.jl

max_stars_repo_name: buiquangbao/hcmus

max_issues_repo_name: buiquangbao/hcmus

max_forks_repo_name: buiquangbao/hcmus


text:
using DifferentialEquations, CSV, DataFrames, DelimitedFiles, StatsBase

function Model!(dy,y,par,t)
  #IFN, ODE 1 parameters
  k11=0 #PR8, RIGI is assumed antagonized
  n=3
  RIGI=1
  k12=par[1]
  k13=par[2]
  k14=par[3]
  #IFN_env, ODE 2 parameters
  k21=par[4]
  tau2=par[5]
  #STATP, ODE 3 parameters
  k31=par[6]
  k32=par[7]
  k33=par[8]
  #IRF7, ODE 4 parameters
  k41=par[9]
  k42=par[10]
  #IRF7P, ODE 5 parameters
  k51=par[11]
  #Cells, ODE 6 parameters
  k61=par[12]
  #Virus, ODE 7 parameters
  k71=par[13]
  k72=par[14]
  k73=par[15]

  #ODE System
  dy[1]=y[6]*(k11*RIGI*y[7]+(k12*y[7]^n)/(k13+y[7]^n)+k14*y[5])-k21*y[1]
  dy[2]=k21*y[1]-tau2*y[2]
  dy[3]=y[6]*(k31*y[2])/(k32+k33*y[2])-0.3*y[3]
  dy[4]=y[6]*(k41*y[3]+k42*y[5])-0.3*y[4]
  dy[5]=y[6]*k51*y[4]-0.3*y[5]
  dy[6]=-k61*y[6]*y[7]
  dy[7]=(k71*y[6]*y[7])/(1+k72*(y[2]*7E-5))-k73*y[7]
end

idx = convert(Matrix,CSV.read("idx.csv",header=false))
chainparams=convert(Matrix,CSV.read("chainparams.csv")) #Read in all sample's parameters
stateNames = ["IFN","IFNe","STATP","IRF7","IRF7P","Live Cells","Virus"]
u0 = [0, 0, 0, 0.72205, 0, 1, 6.9e-8] #Shifted Initial Conditions
tspan = (0,36.0) #Time (start, end)
shift= [7.939, 0, 12.2075, 13.4271, 0, 0, 0] #Noise shift on data
alg = Vern7()   #ODE solver
logSpace=[1,3,4] #Which outputs should be shown in log2 space?
#Listen. This looks like trash, ok? But it takes a smooth time array, your desired time points, and makes a sorted, non-duplicated array
tgraph = rle(sort(vcat(range(tspan[1],tspan[2],length=100),[0.25,0.5,1,1.5,2,4,6,8,12,18,24])))[1]
#tgraph = [0.25,0.5,1,1.5,2,3,4,6,8,9,12,18,24]

function ProbMaker(model,u0,tspan,p)
  prob = ODEProblem(model,u0,tspan,p)
  return prob
end

solMatrix = fill(0.0,(length(idx)),7,length(tgraph))

for i in 1:length(idx)
  p=chainparams[idx[i],:]
  sol = solve(ProbMaker(Model!,u0,tspan,p),alg)
  solMatrix[(i),:,:]=sol(tgraph)[:,:]
  #EXported PARAMeters, LONG/SHORT/BEST version
  CSV.write("exparamtest.csv",DataFrame(solMatrix[(i),:,:]),append=true)
end

#Call python. Pass solution, add chains to Tesla Plot

max_stars_repo_path: Julia Monte Carlo Fitting/TeslaPlot.jl

max_issues_repo_path: Julia Monte Carlo Fitting/TeslaPlot.jl

max_forks_repo_path: Julia Monte Carlo Fitting/TeslaPlot.jl


text:
# ==============================================================================
# Algorithme de Kung (extrait S_N d'un ensemble statique de points S de IR^2)

function getNonDominatedPoints(XFeas, YFeas)
    S = []
    for i=1:length(XFeas)
        push!(S, (XFeas[i] , YFeas[i]) )
    end
    sort!(S, by = x -> x[1])
    SN=[] ; push!(SN, S[1]) ; minYFeas = S[1][2]
    for i=2:length(XFeas)
        if S[i][2] < minYFeas
            push!(SN, S[i]) ; minYFeas = S[i][2]
        end
    end
    return SN
end

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
abstract type Machine <: GeneratorComponent end

@def machine_ports begin
    state_input = [:δ, :ω, :Vf]
    inner_input = [VI_gen_var, VR_gen_var]
end

"""
Parameters of 0-states synchronous machine: Classical Model

# Conmutable structor
```julia
BaseMachine(R, Xd_p, eq_p, MVABase)
```

# Arguments
* `R`::Float64 : Resistance after EMF in machine per unit
* `Xd_p`::Float64 : Reactance after EMF in machine per unit
* `eq_p`::Float64 : Fixed EMF behind the impedance
* `MVABase`::Float64 : Nominal Capacity in MVA

"""
mutable struct BaseMachine <: Machine
  R::Float64
  Xd_p::Float64
  eq_p::Float64
  MVABase::Float64
  n_states::Int64
  states::Vector{Symbol}
  ports::Ports

    function BaseMachine(R::Float64,
                             Xd_p::Float64,
                             eq_p::Float64,
                             MVABase::Float64)
            n_states = 0
            states = Vector{Symbol}()
            @machine_ports

        new(R,
            Xd_p,
            eq_p,
            MVABase,
            n_states,
            states,
            Ports(state_input, inner_input))
    end
end

"""
Parameters of 2-states synchronous machine: One d- and One q-Axis Model
# Conmutable structor
```julia
OneDOneQMachine(R, Xd, Xq, Xd_p, Xq_pp, Td0_p, Tq0_pp, MVABase)
```

# Arguments
* `R`::Float64
* `Xd`::Float64
* `Xq`::Float64
* `Xd_p`::Float64 : Reactance after EMF in machine per unit
* `Xq_pp`::Float64
* `Td0_p`::Float64 : Time constant of transient d-axis voltage
* `Tq0_pp`::Float64 : Time constant of subtransient q-axis voltage
* `MVABase`::Float64 : Nominal Capacity in MVA

"""
mutable struct OneDOneQMachine <: Machine
  R::Float64
  Xd::Float64
  Xq::Float64
  Xd_p::Float64
  Xq_p::Float64
  Td0_p::Float64
  Tq0_p::Float64
  MVABase::Float64
  n_states::Int64
  states::Vector{Symbol}
  ports::Ports

      function OneDOneQMachine(R::Float64,
                                Xd::Float64,
                                Xq::Float64,
                                Xd_p::Float64,
                                Xq_p::Float64,
                                Td0_p::Float64,
                                Tq0_p::Float64,
                                MVABase::Float64)
            n_states = 2
            states = [:eq_p, :ed_p]
            @machine_ports

        new(R,
            Xd,
            Xq,
            Xd_p,
            Xq_p,
            Td0_p,
            Tq0_p,
            MVABase,
            n_states,
            states,
            Ports(state_input, inner_input))
    end
end

"""
Parameters of 6-states synchronous machine: Marconato model
# Conmutable structor
```julia
MarconatoMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, MVABase)
```

# Arguments
* `R`::Float64 : Resistance after EMF in machine per unit
* `Xd`::Float64 : Reactance after EMF in d-axis per unit
* `Xq`::Float64 : Reactance after EMF in q-axis per unit
* `Xd_p`::Float64 : Transient reactance after EMF in d-axis per unit
* `Xq_p`::Float64 : Transient reactance after EMF in q-axis per unit
* `Xd_pp`::Float64 : Subtransient reactance after EMF in d-axis per unit
* `Xq_pp`::Float64 : Subtransient reactance after EMF in q-axis per unit
* `Td0_p`::Float64 : Time constant of transient d-axis voltage
* `Tq0_p`::Float64 : Time constant of transient q-axis voltage
* `Td0_pp`::Float64 : Time constant of subtransient d-axis voltage
* `Tq0_pp`::Float64 : Time constant of subtransient q-axis voltage
* `T_AA`::Float64 : Time constant of d-axis additional leakage
* `MVABase`::Float64 : Nominal Capacity in MVA

"""
mutable struct MarconatoMachine <: Machine
  R::Float64
  Xd::Float64
  Xq::Float64
  Xd_p::Float64
  Xq_p::Float64
  Xd_pp::Float64
  Xq_pp::Float64
  Td0_p::Float64
  Tq0_p::Float64
  Td0_pp::Float64
  Tq0_pp::Float64
  T_AA::Float64
  γd::Float64
  γq::Float64
  MVABase::Float64
  n_states::Int64
  states::Vector{Symbol}
  ports::Ports

      function MarconatoMachine(R::Float64,
                                Xd::Float64,
                                Xq::Float64,
                                Xd_p::Float64,
                                Xq_p::Float64,
                                Xd_pp::Float64,
                                Xq_pp::Float64,
                                Td0_p::Float64,
                                Tq0_p::Float64,
                                Td0_pp::Float64,
                                Tq0_pp::Float64,
                                T_AA::Float64,
                                MVABase::Float64)
            n_states = 6
            states = [:ψq, :ψd, :eq_p, :ed_p, :eq_pp, :ed_pp]
            γd = ( (Td0_pp*Xd_pp)/(Td0_p*Xd_p) )*(Xd-Xd_p)
            γq = ( (Tq0_pp*Xq_pp)/(Tq0_p*Xq_p) )*(Xq-Xq_p)
            @machine_ports

        new(R,
            Xd,
            Xq,
            Xd_p,
            Xq_p,
            Xd_pp,
            Xq_pp,
            Td0_p,
            Tq0_p,
            Td0_pp,
            Tq0_pp,
            T_AA,
            γd,
            γq,
            MVABase,
            n_states,
            states,
            Ports(state_input, inner_input))
    end
end

"""
Parameters of 4-states synchronous machine: Simplified Marconato model
The derivative of stator fluxes (ψd and ψq) is neglected and ωψd = ψd and
ωψq = ψq is assumed (i.e. ω=1.0). This is standard when transmission network
dynamics is neglected.
# Conmutable structor
```julia
SimpleMarconatoMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, MVABase)
```

# Arguments
* `R`::Float64 : Resistance after EMF in machine per unit
* `Xd`::Float64 : Reactance after EMF in d-axis per unit
* `Xq`::Float64 : Reactance after EMF in q-axis per unit
* `Xd_p`::Float64 : Transient reactance after EMF in d-axis per unit
* `Xq_p`::Float64 : Transient reactance after EMF in q-axis per unit
* `Xd_pp`::Float64 : Subtransient reactance after EMF in d-axis per unit
* `Xq_pp`::Float64 : Subtransient reactance after EMF in q-axis per unit
* `Td0_p`::Float64 : Time constant of transient d-axis voltage
* `Tq0_p`::Float64 : Time constant of transient q-axis voltage
* `Td0_pp`::Float64 : Time constant of subtransient d-axis voltage
* `Tq0_pp`::Float64 : Time constant of subtransient q-axis voltage
* `T_AA`::Float64 : Time constant of d-axis additional leakage
* `MVABase`::Float64 : Nominal Capacity in MVA

"""
mutable struct SimpleMarconatoMachine <: Machine
  R::Float64
  Xd::Float64
  Xq::Float64
  Xd_p::Float64
  Xq_p::Float64
  Xd_pp::Float64
  Xq_pp::Float64
  Td0_p::Float64
  Tq0_p::Float64
  Td0_pp::Float64
  Tq0_pp::Float64
  T_AA::Float64
  γd::Float64
  γq::Float64
  MVABase::Float64
  n_states::Int64
  states::Vector{Symbol}
  ports::Ports

      function SimpleMarconatoMachine(R::Float64,
                                      Xd::Float64,
                                      Xq::Float64,
                                      Xd_p::Float64,
                                      Xq_p::Float64,
                                      Xd_pp::Float64,
                                      Xq_pp::Float64,
                                      Td0_p::Float64,
                                      Tq0_p::Float64,
                                      Td0_pp::Float64,
                                      Tq0_pp::Float64,
                                      T_AA::Float64,
                                      MVABase::Float64)
            n_states = 4
            states = [:eq_p, :ed_p, :eq_pp, :ed_pp]
            γd = ( (Td0_pp*Xd_pp)/(Td0_p*Xd_p) )*(Xd-Xd_p)
            γq = ( (Tq0_pp*Xq_pp)/(Tq0_p*Xq_p) )*(Xq-Xq_p)
            @machine_ports

        new(R,
            Xd,
            Xq,
            Xd_p,
            Xq_p,
            Xd_pp,
            Xq_pp,
            Td0_p,
            Tq0_p,
            Td0_pp,
            Tq0_pp,
            T_AA,
            γd,
            γq,
            MVABase,
            n_states,
            states,
            Ports(state_input, inner_input))
    end
end

"""
Parameters of 6-states synchronous machine: Anderson-Fouad model
# Conmutable structor
```julia
AndersonFouadMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, MVABase)
```

# Arguments
* `R`::Float64 : Resistance after EMF in machine per unit
* `Xd`::Float64 : Reactance after EMF in d-axis per unit
* `Xq`::Float64 : Reactance after EMF in q-axis per unit
* `Xd_p`::Float64 : Transient reactance after EMF in d-axis per unit
* `Xq_p`::Float64 : Transient reactance after EMF in q-axis per unit
* `Xd_pp`::Float64 : Subtransient reactance after EMF in d-axis per unit
* `Xq_pp`::Float64 : Subtransient reactance after EMF in q-axis per unit
* `Td0_p`::Float64 : Time constant of transient d-axis voltage
* `Tq0_p`::Float64 : Time constant of transient q-axis voltage
* `Td0_pp`::Float64 : Time constant of subtransient d-axis voltage
* `Tq0_pp`::Float64 : Time constant of subtransient q-axis voltage
* `MVABase`::Float64 : Nominal Capacity in MVA

"""
mutable struct AndersonFouadMachine <: Machine
  R::Float64
  Xd::Float64
  Xq::Float64
  Xd_p::Float64
  Xq_p::Float64
  Xd_pp::Float64
  Xq_pp::Float64
  Td0_p::Float64
  Tq0_p::Float64
  Td0_pp::Float64
  Tq0_pp::Float64
  MVABase::Float64
  n_states::Int64
  states::Vector{Symbol}
  ports::Ports

      function AndersonFouadMachine(R::Float64,
                                Xd::Float64,
                                Xq::Float64,
                                Xd_p::Float64,
                                Xq_p::Float64,
                                Xd_pp::Float64,
                                Xq_pp::Float64,
                                Td0_p::Float64,
                                Tq0_p::Float64,
                                Td0_pp::Float64,
                                Tq0_pp::Float64,
                                MVABase::Float64)
            n_states = 6
            states = [:ψq, :ψd, :eq_p, :ed_p, :eq_pp, :ed_pp]
            @machine_ports

        new(R,
            Xd,
            Xq,
            Xd_p,
            Xq_p,
            Xd_pp,
            Xq_pp,
            Td0_p,
            Tq0_p,
            Td0_pp,
            Tq0_pp,
            MVABase,
            n_states,
            states,
            Ports(state_input, inner_input))
    end
end

"""

Parameters of 4-states simplified Anderson-Fouad (SimpleAFMachine) model.
The derivative of stator fluxes (ψd and ψq) is neglected and ωψd = ψd and
ωψq = ψq is assumed (i.e. ω=1.0). This is standard when transmission network
dynamics is neglected.
If transmission dynamics is considered use the full order Anderson Fouad model.

# Conmutable structor
```julia
SimpleAFMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, MVABase)
```

# Arguments
* `R`::Float64 : Resistance after EMF in machine per unit
* `Xd`::Float64 : Reactance after EMF in d-axis per unit
* `Xq`::Float64 : Reactance after EMF in q-axis per unit
* `Xd_p`::Float64 : Transient reactance after EMF in d-axis per unit
* `Xq_p`::Float64 : Transient reactance after EMF in q-axis per unit
* `Xd_pp`::Float64 : Subtransient reactance after EMF in d-axis per unit
* `Xq_pp`::Float64 : Subtransient reactance after EMF in q-axis per unit
* `Td0_p`::Float64 : Time constant of transient d-axis voltage
* `Tq0_p`::Float64 : Time constant of transient q-axis voltage
* `Td0_pp`::Float64 : Time constant of subtransient d-axis voltage
* `Tq0_pp`::Float64 : Time constant of subtransient q-axis voltage
* `MVABase`::Float64 : Nominal Capacity in MVA

"""
mutable struct SimpleAFMachine <: Machine
  R::Float64
  Xd::Float64
  Xq::Float64
  Xd_p::Float64
  Xq_p::Float64
  Xd_pp::Float64
  Xq_pp::Float64
  Td0_p::Float64
  Tq0_p::Float64
  Td0_pp::Float64
  Tq0_pp::Float64
  MVABase::Float64
  n_states::Int64
  states::Vector{Symbol}
  ports::Ports

      function SimpleAFMachine(R::Float64,
                                Xd::Float64,
                                Xq::Float64,
                                Xd_p::Float64,
                                Xq_p::Float64,
                                Xd_pp::Float64,
                                Xq_pp::Float64,
                                Td0_p::Float64,
                                Tq0_p::Float64,
                                Td0_pp::Float64,
                                Tq0_pp::Float64,
                                MVABase::Float64)
            n_states = 4
            states = [:eq_p, :ed_p, :eq_pp, :ed_pp]
            @machine_ports

        new(R,
            Xd,
            Xq,
            Xd_p,
            Xq_p,
            Xd_pp,
            Xq_pp,
            Td0_p,
            Tq0_p,
            Td0_pp,
            Tq0_pp,
            MVABase,
            n_states,
            states,
            Ports(state_input, inner_input))
    end
end

"""

Parameter of a full order flux stator-rotor model without zero sequence flux in the stator.
The derivative of stator fluxes (ψd and ψq) is NOT neglected. Only one q-axis damping circuit
is considered. All parameters are in machine per unit.

Refer to Chapter 3 of Power System Stability and Control by P. Kundur or Chapter 11 of Power
System Dynamics: Stability and Control, by J. Machowski, J. Bialek and J. Bumby, for more details.
Note that the models are somewhat different (but equivalent) due to the different Park
Transformation used in both books.

# Conmutable structor
```julia
FullMachine(R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_f, L_1d, L_1q)
```

# Arguments
* `R`::Float64 : Stator resistance after EMF in per unit
* `R_f`::Float64 : Field rotor winding resistance in per unit
* `R_1d`::Float64 : Damping rotor winding resistance on d-axis in per unit.
                  This value is denoted as RD in Machowski.
* `R_1q`::Float64 : Damping rotor winding resistance on q-axis in per unit.
                  This value is denoted as RQ in Machowski.
* `L_d`::Float64 : Inductance of fictitious damping that represent the effect
                  of the three-phase stator winding in the d-axis of the rotor, in per unit.
                  This value is denoted as L_ad + L_l in Kundur (and Ld in Machowski).
* `L_q`::Float64 : Inductance of fictitious damping that represent the effect
                  of the three-phase stator winding in the q-axis of the rotor, in per unit.
                  This value is denoted as L_aq + L_l in Kundur.
* `L_ad`::Float64 : Mutual inductance between stator winding and rotor field (and damping)
                    winding inductance on d-axis, in per unit
* `L_aq`::Float64 : Mutual inductance between stator winding and rotor damping winding
                    inductance on q-axis, in per unit
* `L_f1d`::Float64 : Mutual inductance between rotor field winding and rotor damping winding
                    inductance on d-axis, in per unit
* `L_ff`::Float64 : Field rotor winding inductance, in per unit
* `L_1d`::Float64 : Inductance of the d-axis rotor damping circuit, in per unit
* `L_1q`::Float64 : Inductance of the q-axis rotor damping circuit, in per unit
* `MVABase`::Float64 : Nominal Capacity in MVA

"""
mutable struct FullMachine <: Machine
  R::Float64
  R_f::Float64
  R_1d::Float64
  R_1q::Float64
  L_d::Float64
  L_q::Float64
  L_ad::Float64
  L_aq::Float64
  L_f1d::Float64
  L_ff::Float64
  L_1d::Float64
  L_1q::Float64
  MVABase::Float64
  inv_d_fluxlink::Array{Float64, 2}
  inv_q_fluxlink::Array{Float64, 2}
  n_states::Int64
  states::Vector{Symbol}
  ports::Ports

      function FullMachine(R::Float64,
                                R_f::Float64,
                                R_1d::Float64,
                                R_1q::Float64,
                                L_d::Float64,
                                L_q::Float64,
                                L_ad::Float64,
                                L_aq::Float64,
                                L_f1d::Float64,
                                L_ff::Float64,
                                L_1d::Float64,
                                L_1q::Float64,
                                MVABase::Float64)

            inv_d_fluxlink = inv( [ [-L_d L_ad L_ad]; #ψd: Eq 3.127 in Kundur
                                    [-L_ad L_ff L_f1d ]; #ψf: Eq 3.130 in Kundur
                                    [-L_ad L_f1d L_1d ] ]) #ψ1d: Eq 3.131 in Kundur
            inv_q_fluxlink = inv( [ [-L_q L_aq]; #ψq: Eq 3.128 in Kundur
                                    [-L_aq L_1q ] ]) #ψ1q: Eq 3.132 in Kundur
            n_states = 5
            states = [:ψd, :ψq, :ψf, :ψ1d, :ψ1q]
            @machine_ports

        new(R,
            R_f,
            R_1d,
            R_1q,
            L_d,
            L_q,
            L_ad,
            L_aq,
            L_f1d,
            L_ff,
            L_1d,
            L_1q,
            MVABase,
            inv_d_fluxlink,
            inv_q_fluxlink,
            n_states,
            states,
            Ports(state_input, inner_input))
    end
end

"""

Parameter of a full order flux stator-rotor model without zero sequence flux in the stator.
The derivative of stator fluxes (ψd and ψq) is neglected. This is standard when
transmission network dynamics is neglected. Only one q-axis damping circuit
is considered. All per unit are in machine per unit.

Refer to Chapter 3 of Power System Stability and Control by P. Kundur or Chapter 11 of Power
System Dynamics: Stability and Control, by J. Machowski, J. Bialek and J. Bumby, for more details.
Note that the models are somewhat different (but equivalent) due to the different Park
Transformation used in both books.

# Conmutable structor
```julia
SimpleFullMachine(R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_f, L_1d, L_1q)
```

#Arguments
* `R`::Float64 : Stator resistance after EMF in per unit
* `R_f`::Float64 : Field rotor winding resistance in per unit
* `R_1d`::Float64 : Damping rotor winding resistance on d-axis in per unit.
                  This value is denoted as RD in Machowski.
* `R_1q`::Float64 : Damping rotor winding resistance on q-axis in per unit.
                  This value is denoted as RQ in Machowski.
* `L_d`::Float64 : Inductance of fictitious damping that represent the effect
                  of the three-phase stator winding in the d-axis of the rotor, in per unit.
                  This value is denoted as L_ad + L_l in Kundur (and Ld in Machowski).
* `L_q`::Float64 : Inductance of fictitious damping that represent the effect
                  of the three-phase stator winding in the q-axis of the rotor, in per unit.
                  This value is denoted as L_aq + L_l in Kundur.
* `L_ad`::Float64 : Mutual inductance between stator winding and rotor field (and damping)
                    winding inductance on d-axis, in per unit
* `L_aq`::Float64 : Mutual inductance between stator winding and rotor damping winding
                    inductance on q-axis, in per unit
* `L_f1d`::Float64 : Mutual inductance between rotor field winding and rotor damping winding
                    inductance on d-axis, in per unit
* `L_ff`::Float64 : Field rotor winding inductance, in per unit
* `L_1d`::Float64 : Inductance of the d-axis rotor damping circuit, in per unit
* `L_1q`::Float64 : Inductance of the q-axis rotor damping circuit, in per unit
* `MVABase`::Float64 : Nominal Capacity in MVA

"""

mutable struct SimpleFullMachine <: Machine
  R::Float64
  R_f::Float64
  R_1d::Float64
  R_1q::Float64
  L_d::Float64
  L_q::Float64
  L_ad::Float64
  L_aq::Float64
  L_f1d::Float64
  L_ff::Float64
  L_1d::Float64
  L_1q::Float64
  MVABase::Float64
  inv_d_fluxlink::Array{Float64, 2}
  inv_q_fluxlink::Array{Float64, 2}
  n_states::Int64
  states::Vector{Symbol}
  ports::Ports

      function SimpleFullMachine(R::Float64,
                                R_f::Float64,
                                R_1d::Float64,
                                R_1q::Float64,
                                L_d::Float64,
                                L_q::Float64,
                                L_ad::Float64,
                                L_aq::Float64,
                                L_f1d::Float64,
                                L_ff::Float64,
                                L_1d::Float64,
                                L_1q::Float64,
                                MVABase::Float64)

            inv_d_fluxlink = inv( [ [-L_d L_ad L_ad]; #ψd: Eq 3.127 in Kundur
                                    [-L_ad L_ff L_f1d ]; #ψf: Eq 3.130 in Kundur
                                    [-L_ad L_f1d L_1d ] ]) #ψ1d: Eq 3.131 in Kundur
            inv_q_fluxlink = inv( [ [-L_q L_aq]; #ψq: Eq 3.128 in Kundur
                                    [-L_aq L_1q ] ]) #ψ1q: Eq 3.132 in Kundur
            n_states = 3
            states = [:ψf, :ψ1d, :ψ1q]
            @machine_ports

        new(R,
            R_f,
            R_1d,
            R_1q,
            L_d,
            L_q,
            L_ad,
            L_aq,
            L_f1d,
            L_ff,
            L_1d,
            L_1q,
            MVABase,
            inv_d_fluxlink,
            inv_q_fluxlink,
            n_states,
            states,
            Ports(state_input, inner_input))
    end
end


text:
N(x) = sum(x)

function fluxes(model)
    terms = map(enumerate(model.Δ)) do (i, δ)
        inn, out = δ
        term = inn
        if term.op == +
            term = prod(inn.args)
        end
        body, args = Petri.funcbody(term)
        t = :(param[$i]*$body)
        deg = length(t.args[3].args)-2
        t = :($t / N(state)^$(deg))
        t = simplify(t)
        outterms = out.op == ( + ) ? out.args : [out]
        changes = map(outterms) do o
            p = o.op == ( * ) ? (o.args[2].op,o.args[1]) : (o.op, 1)
            var, coeff = p[1], coeffvalue(p[2])
            var=>:($coeff * $t)
        end
        innterms = inn.op == ( + ) ? inn.args : [inn]
        decreases = map(innterms) do o
            p = o.op == ( * ) ? (o.args[2].op,o.args[1]) : (o.op, 1)
            var, coeff = p[1], coeffvalue(p[2])
            var=>:((-1/$coeff) * $t)
        end
        return union(changes, decreases)
    end |> flatten |> s->collect(Dict(), s)
    [:(du.$(Symbol(k)) = +($(v...))) for (k,v) in terms]
end

odefunc(m, prefix::Symbol) = funckit(gensym(prefix),
                                     (:du, :state, :param, :t),
                                     quotesplat(fluxes(m)))

function mk_function(m::Model)
    fex = odefunc(m, :f)
    arglist = fex.args[2].args[1].args[2:end]
    body = fex.args[2].args[2].args[2]
    f = mk_function(@__MODULE__, arglist, [], body)
    g = (du, u, p, t) -> f(du,  u, p, t)
    return g
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayHttpListenerPropertiesFormat <: SwaggerModel
    frontendIPConfiguration::Any # spec type: Union{ Nothing, SubResource } # spec name: frontendIPConfiguration
    frontendPort::Any # spec type: Union{ Nothing, SubResource } # spec name: frontendPort
    protocol::Any # spec type: Union{ Nothing, ApplicationGatewayProtocol } # spec name: protocol
    hostName::Any # spec type: Union{ Nothing, String } # spec name: hostName
    sslCertificate::Any # spec type: Union{ Nothing, SubResource } # spec name: sslCertificate
    requireServerNameIndication::Any # spec type: Union{ Nothing, Bool } # spec name: requireServerNameIndication
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    customErrorConfigurations::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayCustomError} } # spec name: customErrorConfigurations
    firewallPolicy::Any # spec type: Union{ Nothing, SubResource } # spec name: firewallPolicy
    hostNames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: hostNames

    function ApplicationGatewayHttpListenerPropertiesFormat(;frontendIPConfiguration=nothing, frontendPort=nothing, protocol=nothing, hostName=nothing, sslCertificate=nothing, requireServerNameIndication=nothing, provisioningState=nothing, customErrorConfigurations=nothing, firewallPolicy=nothing, hostNames=nothing)
        o = new()
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        setfield!(o, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("frontendPort"), frontendPort)
        setfield!(o, Symbol("frontendPort"), frontendPort)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("hostName"), hostName)
        setfield!(o, Symbol("hostName"), hostName)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("sslCertificate"), sslCertificate)
        setfield!(o, Symbol("sslCertificate"), sslCertificate)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("requireServerNameIndication"), requireServerNameIndication)
        setfield!(o, Symbol("requireServerNameIndication"), requireServerNameIndication)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("customErrorConfigurations"), customErrorConfigurations)
        setfield!(o, Symbol("customErrorConfigurations"), customErrorConfigurations)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("firewallPolicy"), firewallPolicy)
        setfield!(o, Symbol("firewallPolicy"), firewallPolicy)
        validate_property(ApplicationGatewayHttpListenerPropertiesFormat, Symbol("hostNames"), hostNames)
        setfield!(o, Symbol("hostNames"), hostNames)
        o
    end
end # type ApplicationGatewayHttpListenerPropertiesFormat

const _property_map_ApplicationGatewayHttpListenerPropertiesFormat = Dict{Symbol,Symbol}(Symbol("frontendIPConfiguration")=>Symbol("frontendIPConfiguration"), Symbol("frontendPort")=>Symbol("frontendPort"), Symbol("protocol")=>Symbol("protocol"), Symbol("hostName")=>Symbol("hostName"), Symbol("sslCertificate")=>Symbol("sslCertificate"), Symbol("requireServerNameIndication")=>Symbol("requireServerNameIndication"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("customErrorConfigurations")=>Symbol("customErrorConfigurations"), Symbol("firewallPolicy")=>Symbol("firewallPolicy"), Symbol("hostNames")=>Symbol("hostNames"))
const _property_types_ApplicationGatewayHttpListenerPropertiesFormat = Dict{Symbol,String}(Symbol("frontendIPConfiguration")=>"SubResource", Symbol("frontendPort")=>"SubResource", Symbol("protocol")=>"ApplicationGatewayProtocol", Symbol("hostName")=>"String", Symbol("sslCertificate")=>"SubResource", Symbol("requireServerNameIndication")=>"Bool", Symbol("provisioningState")=>"ProvisioningState", Symbol("customErrorConfigurations")=>"Vector{ApplicationGatewayCustomError}", Symbol("firewallPolicy")=>"SubResource", Symbol("hostNames")=>"Vector{String}")
Base.propertynames(::Type{ ApplicationGatewayHttpListenerPropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayHttpListenerPropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayHttpListenerPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayHttpListenerPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayHttpListenerPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayHttpListenerPropertiesFormat[property_name]

function check_required(o::ApplicationGatewayHttpListenerPropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayHttpListenerPropertiesFormat }, name::Symbol, val)
end


text:
module mfes1
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim, RefShapeInterval
using Elfel.FElements: FE, refshape, FEH1_L2, FEH1_T3, FEH1_Q4
using Elfel.FElements: bfun, bfungradpar, nfeatofdim, ndofperfeat
using Test
function test()
    fe = FEH1_T3()
    @test ndofperfeat(fe, 0) == 1
    @test ndofperfeat(fe, 1) == 0
    @test ndofperfeat(fe, 2) == 0
    @test ndofperfeat(fe, 3) == 0

    @test refshape(fe) == RefShapeTriangle
    @test manifdim(refshape(fe)) == 2
    @test nfeatofdim(fe, 0) == 3
    @test nfeatofdim(fe, 1) == 3
    @test nfeatofdim(fe, 2) == 1
    @test nfeatofdim(fe, 3) == 0
    @test isapprox(bfun(fe, [1/3, 1/3]), [0.3333333333333334, 0.3333333333333333, 0.3333333333333333])
    g = bfungradpar(fe, [1/3, 1/3])
    @test isapprox(g[1], [-1.; -1.]')
    @test isapprox(g[2], [+1.;  0.]')
    @test isapprox(g[3], [0.; +1.]')

    e = FEH1_T3()
    @test ndofperfeat(fe, 0) == true
    @test ndofperfeat(fe, 1) == 0
    @test ndofperfeat(fe, 2) == 0
    @test ndofperfeat(fe, 3) == 0
    @test refshape(fe) == RefShapeTriangle
    @test manifdim(refshape(fe)) == 2
    @test isapprox(bfun(fe, [1/3, 1/3]), [0.3333333333333334, 0.3333333333333333, 0.3333333333333333])
    g = bfungradpar(fe, [1/3, 1/3])
    @test isapprox(g[1], [-1.; -1.]')
    @test isapprox(g[2], [+1.;  0.]')
    @test isapprox(g[3], [0.; +1.]')

    fe = FEH1_L2()
    @test ndofperfeat(fe, 0) == true
    @test ndofperfeat(fe, 1) == 0
    @test ndofperfeat(fe, 2) == 0
    @test ndofperfeat(fe, 3) == 0
    @test refshape(fe) == RefShapeInterval
    @test manifdim(refshape(fe)) == 1
    @test isapprox(bfun(fe, [1/3]), [0.3333333333333334, 2*0.3333333333333333])
    g = bfungradpar(fe, [1/3])
    @test isapprox(g[1], [-1.0/2])
    @test isapprox(g[2], [+1.0/2])
end
end
using .mfes1
mfes1.test()

module mfes2
using Elfel
using Elfel.RefShapes: manifdim, RefShapeTetrahedron
using Elfel.FElements: FE, refshape, FEH1_T4
using Elfel.FElements: bfun, bfungradpar, nfeatofdim, ndofperfeat
using Test
function test()
    fe = FEH1_T4()
    @test ndofperfeat(fe, 0) == 1
    @test ndofperfeat(fe, 1) == 0
    @test ndofperfeat(fe, 2) == 0
    @test ndofperfeat(fe, 3) == 0

    @test refshape(fe) == RefShapeTetrahedron
    @test manifdim(refshape(fe)) == 3
    @test nfeatofdim(fe, 0) == 4
    @test nfeatofdim(fe, 1) == 6
    @test nfeatofdim(fe, 2) == 4
    @test nfeatofdim(fe, 3) == 1
    @test isapprox(bfun(fe, [1/4, 1/4, 1/4]), [1/4, 1/4, 1/4, 1/4])
    g = bfungradpar(fe, [1/4, 1/4, 1/4])
    @test isapprox(g[1], [-1. -1. -1.])
    @test isapprox(g[2], [+1.  0. 0.])
    @test isapprox(g[3], [0. +1.  0.])
    @test isapprox(g[4], [0. 0. +1. ])
true
end
end
using .mfes2
mfes2.test()

module mfes3
using Elfel
using Elfel.RefShapes: manifdim, RefShapeSquare
using Elfel.FElements: FE, refshape, FEL2_Q4
using Elfel.FElements: bfun, bfungradpar, nfeatofdim, ndofperfeat
using Test
function test()
    fe = FEL2_Q4()
    @test ndofperfeat(fe, 0) == 0
    @test ndofperfeat(fe, 1) == 0
    @test ndofperfeat(fe, 2) == 1
    @test ndofperfeat(fe, 3) == 0

    @test refshape(fe) == RefShapeSquare
    @test manifdim(refshape(fe)) == 2
    @test nfeatofdim(fe, 0) == 4
    @test nfeatofdim(fe, 1) == 4
    @test nfeatofdim(fe, 2) == 1
    @test nfeatofdim(fe, 3) == 0
    @test isapprox(bfun(fe, [1/4, 1/4]), [1.0])
    g = bfungradpar(fe, [1/4, 1/4])
    @test isapprox(g[1], [0. 0. ])
true
end
end
using .mfes3
mfes3.test()

module mfes4
using Elfel
using Elfel.RefShapes: manifdim, RefShapeTetrahedron, RefShapeTriangle
using Elfel.FElements: FE, refshape, FEL2_T4, FEL2_T3
using Elfel.FElements: bfun, bfungradpar, nfeatofdim, ndofperfeat
using Test
function test()
    fe = FEL2_T4()
    @test ndofperfeat(fe, 0) == 0
    @test ndofperfeat(fe, 1) == 0
    @test ndofperfeat(fe, 2) == 0
    @test ndofperfeat(fe, 3) == 1

    @test refshape(fe) == RefShapeTetrahedron
    @test manifdim(refshape(fe)) == 3
    @test nfeatofdim(fe, 0) == 4
    @test nfeatofdim(fe, 1) == 6
    @test nfeatofdim(fe, 2) == 4
    @test nfeatofdim(fe, 3) == 1
    @test isapprox(bfun(fe, [1/4, 1/4, 1/4]), [1.0])
    g = bfungradpar(fe, [1/4, 1/4, 1/4])
    @test isapprox(g[1], [0.0 0.0 0.0])

    fe = FEL2_T3()
    @test ndofperfeat(fe, 0) == 0
    @test ndofperfeat(fe, 1) == 0
    @test ndofperfeat(fe, 2) == 1
    @test ndofperfeat(fe, 3) == 0

    @test refshape(fe) == RefShapeTriangle
    @test manifdim(refshape(fe)) == 2
    @test nfeatofdim(fe, 0) == 3
    @test nfeatofdim(fe, 1) == 3
    @test nfeatofdim(fe, 2) == 1
    @test nfeatofdim(fe, 3) == 0
    @test isapprox(bfun(fe, [1/4, 1/4]), [1.0])
    g = bfungradpar(fe, [1/4, 1/4])
    @test isapprox(g[1], [0.0 0.0])
true
end
end
using .mfes4
mfes4.test()


text:
###############################################################################
# Reset and Clear
###############################################################################
export @clearsnaps, @clearallsnaps, @resetallsnaps

brandnewsnaps!() = begin
    global _num_trace_locations = 0
    empty!(happysnaps)
    empty!(location_counts)
    empty!(watched_exprs)
end

macro resetallsnaps()
    :(RickTracy.brandnewsnaps!()) |> esc
end

macro clearallsnaps()
    :(empty!(RickTracy.happysnaps)) |> esc
end

clearsnaps(exprstr) = begin
    #find snaps that match key, and remove them from the happysnaps vector
    filter!((st)->st.exprstr == exprstr, happysnaps) #slow
end

macro clearsnaps(exprs...)
    res = :()
    for expr in exprs
        exprstr = "$expr" #expr as a string
        res = :($res; RickTracy.clearsnaps(exprstr))
    end
    res |> esc
end

macro clearunwatch(exprs...)
    :(@unwatch exprs; @clearsnaps(exprs)) |> esc
end

max_stars_repo_name: JobJob/RickTracy


text:
using SIUnits: SIQuantity, SIUnit, unit
using SIUnits.ShortUnits
using LatBo

include("SIUnitsUnits.jl")
using LatBoTests.Units: as, from

facts("Physical to and from LB units using SIUnits") do
    is_dimensionless{T<:Number}(::Type{T}) = !(T <: SIQuantity) && !(T <: SIUnit)
    is_dimensionless(x::Number) = is_dimensionless(typeof(x))
    is_dimensionless(x::Array) = is_dimensionless(eltype(x))

    δt, δx, δm = 0.1, 2., 30.
    siunits = LatBoTests.Units.LBUnits{Float64}(δt * s, δx * m, δm * kg)

    context("P to LB") do
        @fact as(siunits, 0.1m) --> is_dimensionless
        @fact as(siunits, 2m) --> roughly(2 / δx)
        @fact as(siunits, 1m) --> roughly(1 / δx)
        @fact as(siunits, 1/m) --> roughly(δx)
        @fact as(siunits, 1s) --> is_dimensionless
        @fact as(siunits, 1s) --> roughly(1 / δt)
        @fact as(siunits, 1m/s) --> roughly(δt/δx)
        @fact as(siunits, 0.5kg/m^3/s) --> roughly(0.5δt*δx^3/δm)
        @fact as(siunits, [0.5, 0.6]*(kg*m^-3/s)) --> roughly([0.5, 0.6]δt*δx^3/δm)
        @fact as(siunits, [0.5, 0.6]*(kg*m^-3/s)) --> is_dimensionless
        @fact_throws as(siunits, 0.1Pascal * mol)
    end

    context("LB to P") do
        @fact from(siunits, m) --> roughly(siunits.δx)
        @fact unit(from(siunits, 2, m)) --> m
        @fact from(siunits, 2, m) --> roughly(2.siunits.δx)
        some_dim = siunits.δm * siunits.δx^-3 * siunits.δt^-1
        @fact from(siunits, 0.5, kg/m^3/s) --> roughly(0.5 * some_dim)
        @fact from(siunits, [0.5, 0.6], kg*m^-3/s) --> roughly([0.5, 0.6] * some_dim)
    end

    context("P to LB to P") do
        for value in [1m, 2m, 3kg/s/m^2]
            @fact from(siunits, as(siunits, value), unit(value)) --> roughly(value)
        end
    end
end

facts("Physical to and from LB units") do
    δt, δx, δm = 0.1, 2., 30.
    lbunits = LatBo.Units.LBUnits{Float64}(δt, δx, δm)
    siunits = LatBoTests.Units.LBUnits{Float64}(δt * s, δx * m, δm * kg)

    equivs = [
        (:time, s), (:length, m), (:weight, kg), (:density, 1kg*m^-3),
        (:velocity, 1m/s), (:pressure, Pa), (:viscosity, 1Pa * s)
    ]
    goaround(dim::Symbol, x) = LatBo.Units.physical(
        lbunits, dim, LatBo.Units.dimensionless(lbunits, dim, x))
    for (lb, si) in equivs
        @fact LatBo.Units.dimensionless(lbunits, lb, 2) --> roughly(as(siunits, 2 * si))
        r = 0.5 * randn(10, 5) + 0.5
        @fact LatBo.Units.dimensionless(lbunits, lb, r) --> roughly(as(siunits, r * si))

        @fact goaround(lb, 2) --> roughly(2)
        @fact goaround(lb, r) --> roughly(r)
    end

    @fact_throws LatBo.Units.dimensionles(lbunits, :notaunit, 2)
end


text:
using GPX
using TimeZones
using LightXML: XMLDocument, save_file

author = GPXAuthor("S. Celles")

metadata = GPXMetadata(
    name="07/11/2019 LFBI (09:32) LFBI (11:34)",
    author=author,
    time=ZonedDateTime("2019-01-01T00:00:0.000+00:00"), # ZonedDateTime("2019-01-01T00:00:00.000Z") when https://github.com/JuliaTime/TimeZones.jl/pull/227 will be merged
)

gpx = GPXDocument(metadata)

track = new_track(gpx)
track_segment = new_track_segment(track)
point = GPXPoint(0.0, 0.0, 0.0, ZonedDateTime("2019-11-09T09:14:17.473+00:00"), "")
push!(track_segment, point)

# println(gpx)

xdoc = XMLDocument(gpx)

println(xdoc)

fname = "generated.gpx"
save_file(xdoc, fname)
println("GPX file saved to \"$fname\"")


text:
using Test
using MathOptInterface
import Clp

const MOI  = MathOptInterface

const OPTIMIZER = Clp.Optimizer()
MOI.set(OPTIMIZER, MOI.Silent(), true)

@testset "SolverName" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "Clp"
end

@testset "supports_default_copy_to" begin
    @test !MOI.Utilities.supports_allocate_load(OPTIMIZER, false)
    @test !MOI.Utilities.supports_allocate_load(OPTIMIZER, true)
    @test !MOI.Utilities.supports_default_copy_to(OPTIMIZER, false)
    @test !MOI.Utilities.supports_default_copy_to(OPTIMIZER, true)
end

const CACHED = MOI.Utilities.CachingOptimizer(
    MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    OPTIMIZER
)

const BRIDGED = MOI.Bridges.full_bridge_optimizer(CACHED, Float64)

const CONFIG = MOI.Test.TestConfig(
    dual_objective_value = false,
)

@testset "basic_constraint_tests" begin
    MOI.Test.basic_constraint_tests(CACHED, CONFIG)
end

@testset "Unit Tests" begin
    MOI.Test.unittest(BRIDGED, CONFIG, [
        # Not supported by upstream.
        "number_threads",

        # Tests that require integer variables
        "solve_integer_edge_cases",
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3",
        "solve_objbound_edge_cases",

        # Tests that require quadratic objective / constraints
        "solve_qcp_edge_cases",
        "solve_qp_edge_cases",
        "delete_soc_variables",
    ])
end

@testset "Linear tests" begin
    MOI.Test.contlineartest(BRIDGED, CONFIG, [
        # The linear12 test requires the InfeasibilityCertificate for variable
        # bounds. These are available through C++, but not via the C interface.
        "linear12",

        # MOI.VariablePrimalStart not supported.
        "partial_start"
    ])
end

@testset "ModelLike" begin
    @testset "nametest" begin
        MOI.Test.nametest(BRIDGED)
    end
    @testset "validtest" begin
        MOI.Test.validtest(BRIDGED)
    end
    @testset "emptytest" begin
        MOI.Test.emptytest(BRIDGED)
    end
end

@testset "Nonexistant unbounded ray" begin
    MOI.empty!(BRIDGED)
    x = MOI.add_variables(BRIDGED, 5)
    MOI.set(
        BRIDGED,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0)
    )
    MOI.set(BRIDGED, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.optimize!(BRIDGED)
    status = MOI.get(BRIDGED, MOI.TerminationStatus())
    @test status == MOI.DUAL_INFEASIBLE
end

@testset "RawParameter" begin
    model = Clp.Optimizer()
    MOI.set(model, MOI.RawParameter("LogLevel"), 1)
    @test MOI.get(model, MOI.RawParameter("LogLevel")) == 1
    @test MOI.get(model, MOI.RawParameter(:LogLevel)) == 1
    MOI.set(model, MOI.RawParameter(:LogLevel), 2)
    @test MOI.get(model, MOI.RawParameter("LogLevel")) == 2
    @test MOI.get(model, MOI.RawParameter(:LogLevel)) == 2

    MOI.set(model, MOI.RawParameter("SolveType"), 1)
    @test MOI.get(model, MOI.RawParameter("SolveType")) == 1
    @test MOI.get(model, MOI.RawParameter(:SolveType)) == 1
    MOI.set(model, MOI.RawParameter("SolveType"), 4)
    @test MOI.get(model, MOI.RawParameter("SolveType")) == 4
    @test MOI.get(model, MOI.RawParameter(:SolveType)) == 4

    MOI.set(model, MOI.RawParameter("PresolveType"), 1)
    @test MOI.get(model, MOI.RawParameter("PresolveType")) == 1
    @test MOI.get(model, MOI.RawParameter(:PresolveType)) == 1
    MOI.set(model, MOI.RawParameter("PresolveType"), 0)
    @test MOI.get(model, MOI.RawParameter("PresolveType")) == 0
    @test MOI.get(model, MOI.RawParameter(:PresolveType)) == 0
end


text:
#=
    This file contains all of the functions that a solver needs to implement
    in order to use LQOI.

       min/max: c'x + x'Qx
    subject to:
                # Variable bounds
                abᵢ <=  xᵢ              <= bbᵢ, i=1,2,...Nx
                # Linear Constraints
                llᵢ <= alᵢ'x             <= ulᵢ, i=1,2,...Nl
                # Quadratic Constraints
                lqᵢ <= aqᵢ' x + x'Qqᵢ' x <= uqᵢ, i=1,2,...Nq
                # SOS1, SOS2 constraints
                # Binary Constraints
                xᵢ ∈ {0, 1}
                # Integer Constraints
                xᵢ ∈ Z
=#

"""
    backend_type(m, ::MOI.AbstractSet)::CChar

An overloadable type for dispatching the appropriate types to the backends.

For example, GLPK.jl uses `'E'` for `a'x=b` constraints, where as Gurobi.jl uses
`==`.

Three are special cases:
 - `Val{:Continuous}`: for the type of a continuous variable
 - `Val{:Upperbound}`: for the upper bound of a variable
 - `Val{:Lowerbound}`: for the lower bound of a variable

### Defaults

    MOI.GreaterThan  - 'G'
    MOI.LessThan     - 'L'
    MOI.EqualTo      - 'E'

    MOI.Zeros        - 'E'
    MOI.Nonpositives - 'L'
    MOI.Nonnegatives - 'G'

    MOI.ZeroOne        - 'B'
    MOI.Integer        - 'I'
    MOI.Semicontinuous - 'S'
    MOI.Semiinteger    - 'N'

    MOI.SOS1 - :SOS1  # '1'
    MOI.SOS2 - :SOS2  # '2'

    Val{:Continuous} - 'C'
    Val{:Upperbound} - 'U'
    Val{:Lowerbound} - 'L'
"""
backend_type(m::LinQuadOptimizer, ::MOI.GreaterThan{T}) where T = Cchar('G')
backend_type(m::LinQuadOptimizer, ::MOI.LessThan{T}) where T    = Cchar('L')
backend_type(m::LinQuadOptimizer, ::MOI.EqualTo{T}) where T     = Cchar('E')
# Implemented separately
# backend_type(m::LinQuadOptimizer, ::MOI.Interval{T}) where T    = Cchar('R')

backend_type(m::LinQuadOptimizer, ::MOI.Zeros)        = Cchar('E')
backend_type(m::LinQuadOptimizer, ::MOI.Nonpositives) = Cchar('L')
backend_type(m::LinQuadOptimizer, ::MOI.Nonnegatives) = Cchar('G')

backend_type(m::LinQuadOptimizer, ::MOI.ZeroOne) = Cchar('B')
backend_type(m::LinQuadOptimizer, ::MOI.Integer) = Cchar('I')

backend_type(m::LinQuadOptimizer, ::MOI.SOS1{T}) where T = :SOS1  # Cchar('1')
backend_type(m::LinQuadOptimizer, ::MOI.SOS2{T}) where T = :SOS2  # Cchar('2')

backend_type(m::LinQuadOptimizer, ::MOI.Semicontinuous{T}) where T = Cchar('S')
backend_type(m::LinQuadOptimizer, ::MOI.Semiinteger{T}) where T    = Cchar('N')

backend_type(m::LinQuadOptimizer, ::Val{:Continuous}) = Cchar('C')
backend_type(m::LinQuadOptimizer, ::Val{:Upperbound}) = Cchar('U')
backend_type(m::LinQuadOptimizer, ::Val{:Lowerbound}) = Cchar('L')

"""
    LinearQuadraticModel(M::Type{<:LinQuadOptimizer}, env)

Initializes a model given a model type `M` and an `env` that might be a `nothing`
for some solvers.
"""
function LinearQuadraticModel end

"""
    supported_constraints(m)::Vector{
        Tuple{MOI.AbstractFunction, MOI.AbstractSet}
    }

Get a list of supported constraint types in the model `m`.

For example, `[(LQOI.Linear, LQOI.EQ)]`
"""
function supported_constraints end

"""
    supported_objectives(m)::Vector{MOI.AbstractScalarFunction}

Get a list of supported objective types in the model `m`.

For example, `[LQOI.Linear, LQOI.Quad]`
"""
function supported_objectives end

# Constraints

"""
    change_variable_bounds!(m, cols::Vector{Int}, values::Vector{Float64}, senses::Vector)

Change the bounds of the variable. The sense of the upperbound
is given by `backend_type(m, Val{:Upperbound}())`. The sense
of the lowerbound is given by `backend_type(m, Val{:Lowerbound}())`
"""
function change_variable_bounds! end

"""
    get_variable_lowerbound(m, col::Int)::Float64

Get the lower bound of the variable in 1-indexed column `col` of the model `m`.
"""
function get_variable_lowerbound end

"""
    get_variable_upperbound(m, col::Int)::Float64

Get the upper bound of the variable in 1-indexed column `col` of the model `m`.
"""
function get_variable_upperbound end

"""
    get_number_linear_constraints(m)::Int

Get the number of linear constraints in the model `m`.
"""
function get_number_linear_constraints end

"""
    add_linear_constraints!(m, A::CSRMatrix{Float64},
        sense::Vector{Cchar}, rhs::Vector{Float64})::Nothing

Adds linear constraints of the form `Ax (sense) rhs` to the model `m`.

`sense` and `rhs` contain one element for each row in `A`.
The `sense` is given by `backend_type(m, set)`.

Ranged constraints (`set=MOI.Interval`) should be added via `add_ranged_constraint!`
instead.

See also: `LinQuadOptInterface.CSRMatrix`.
"""
function add_linear_constraints! end

"""
    add_ranged_constraints!(m, A::CSRMatrix{Float64},
        lowerbound::Vector{Float64}, upperbound::Vector{Float64})

Adds linear constraints of the form `lowerbound <= Ax <= upperbound` to the
model `m`.

This is a special case compared to standard `add_linear_constraints!` since it
is often implemented via multiple API calls.
"""
function add_ranged_constraints! end

"""
    modify_ranged_constraints!(m, rows::Vector{Int}, lowerbound::Vector{Float64}, upperbound::Vector{Float64})

Modify the lower and upperbounds of a ranged constraint in the model `m`.

This is a special case compared to standard the `change_rhs_coefficient!` since it
is often implemented via multiple API calls.
"""
function modify_ranged_constraints! end

"""
    get_rhs(m, row::Int)::Float64

Get the right-hand side of the linear constraint in the 1-indexed row `row` in
the model `m`.
"""
function get_rhs end

"""
    get_range(m, row::Int)::Tuple{Float64,Float64}

Get the range which the constraint `row` belongs to. The output of the function is the tuple `lowerbound, upperbound` of bounds: `lowerbound <= a'x < = upperbound`
"""
function get_range end

"""
    get_linear_constraint(m, row::Int)::Tuple{Vector{Int}, Vector{Float64}}

Get the linear component of the constraint in the 1-indexed row `row` in
the model `m`. Returns a tuple of `(cols, vals)`.
"""
function get_linear_constraint end

"""
    change_matrix_coefficient!(m, row, col, coef)

Set the linear coefficient of the variable in column `col`, constraint `row` to
`coef`.
"""
function change_matrix_coefficient! end

"""
    change_matrix_coefficients!(m, rows, cols, coefs)

Change multiple linear coefficients simultaneously. Sets the linear coefficient
of the variable at column `cols[i]` in constraint `rows[i]` to `coefs[i]` for
`i in 1:length(rows)`. Requires that `length(rows) == length(cols) == length(coefs)`.

By default, this method just calls
`change_matrix_coefficient!(m, rows[i], cols[i], coefs[i])` repeatedly, but
some solver interfaces may offer more efficient implementations.
"""
function change_matrix_coefficients!(m, rows, cols, coefs)
    @boundscheck((axes(rows) == axes(cols) == axes(coefs)) || throw(DimensionMismatch()))
    for i in eachindex(rows)
        change_matrix_coefficient!(m, rows[i], cols[i], coefs[i])
    end
end

"""
change_objective_coefficient!(m, col, coef)

Set the linear coefficient of the variable in column `col` to `coef` in the objective function.
"""
function change_objective_coefficient! end

"""
change_rhs_coefficient!(m, row, coef)

Set the rhs of the constraint in row `row` to `coef`.
"""
function change_rhs_coefficient! end

"""
    delete_linear_constraints!(m, start_row::Int, end_row::Int)::Nothing

Delete the linear constraints `start_row`, `start_row+1`, ..., `end_row` from
the model `m`.
"""
function delete_linear_constraints! end

"""
    delete_quadratic_constraints!(m, start_row::Int, end_row::Int)::Nothing

Delete the quadratic constraints `start_row`, `start_row+1`, ..., `end_row` from
the model `m`.
"""
function delete_quadratic_constraints! end

"""
    change_variable_types(m, cols::Vector{Int}, types):Nothing

Change the variable types. Type is the output of one of:
 - `backend_type(m, ::ZeroOne)`, for binary variables;
 - `backend_type(m, ::Integer)`, for integer variables; and
 - `backend_type(m, Val{:Continuous}())`, for continuous variables.
"""
function change_variable_types! end

"""
    change_linear_constraint_sense!(m, rows::Vector{Int}, sense::Vector{Symbol})::Nothing

Change the sense of the linear constraints in `rows` to `sense`.

`sense` is the output of `backend_type(m, set)`, where `set`
is the corresponding set for the row `rows[i]`.

`Interval` constraints require a call to `change_range_value!`.
"""
function change_linear_constraint_sense! end

"""
    make_problem_type_integer(m)::Nothing

If an explicit call is needed to change the problem type integer (e.g., CPLEX).
"""
function make_problem_type_integer(m::LinQuadOptimizer)
    nothing  # default
end

"""
    make_problem_type_continuous(m)::Nothing

If an explicit call is needed to change the problem type continuous (e.g., CPLEX).
"""
function make_problem_type_continuous(m::LinQuadOptimizer)
    nothing  # default
end

"""
    add_sos_constraint!(m, cols::Vector{Int}, vals::Vector{Float64}, typ::Symbol)::Nothing

Add the SOS constraint to the model `m`. `typ` is either `:SOS1` or `:SOS2`.
"""
function add_sos_constraint! end

"""
    delete_sos!(m, start_idx::Int, end_idx::Int)::Nothing

Delete the SOS constraints `start_idx`, `start_idx+1`, ..., `end_idx` from
the model `m`.
"""
function delete_sos! end

"""
    get_sos_constraint(m, idx::Int)::Tuple{Vector{Int}, Vector{Float64}, Symbol}

Get the SOS constraint `idx` from the model `m`. Returns the triplet
    `(cols, vals, typ)`.
"""
function get_sos_constraint end

"""
    get_number_quadratic_constraints(m)::Int

Get the number of quadratic constraints in the model `m`.
"""
function get_number_quadratic_constraints end

"""
    add_quadratic_constraint!(m, cols::Vector{Int}, coefs::Vector{Float64}, rhs::Float64,
        sense, I::Vector{Int}, J::Vector{Int}, V::Vector{Float64})::Nothing

Add a quadratic constraint `a'x + 0.5 x' Q x`.
See `add_linear_constraints!` for information of linear component.
Arguments `(I,J,V)` given in triplet form for the Q matrix in `0.5 x' Q x`.
"""
function add_quadratic_constraint! end

"""
    set_constant_objective!(m, value)::Nothing

Set the constant (i.e. offset) component of the objective function to the given
value.

Solver interfaces that overload this behavior (e.g. to pass that constant
objective to the solver itself) must also overload `get_constant_objective(m)`.
"""
function set_constant_objective!(m::LinQuadOptimizer, value)
    m.objective_constant = value
    return
end

"""
    set_quadratic_objective!(m, I::Vector{Int}, J::Vector{Int}, V::Vector{Float64})::Nothing

Set the quadratic component of the objective. Arguments given in triplet form
for the Q matrix in `0.5 x' Q x`.
"""
function set_quadratic_objective! end

"""
    get_quadratic_constraint(m, row::Int)::Tuple{Vector{Int}, Vector{Float64}, SparseMatrixCSC{Float64,Int64}}

Get the linear and quadratic components of the constraint in the 1-indexed row `row` in
the model `m`. Returns a tuple of `(lin_cols, lin_vals, Q)`.
Where `Q` represents the matrix in CSC format.
"""
function get_quadratic_constraint end

"""
    get_quadratic_rhs(m, row::Int)::Float64

Get the right hand-side term of quadratic constraint in row `row` in model `m`.
"""
function get_quadratic_rhs end

"""
    set_linear_objective!(m, cols::Vector{Int}, coefs::Vector{Float64})::Nothing

Set the linear component of the objective.
"""
function set_linear_objective! end

"""
    change_objective_sense!(m, sense::Symbol)::Nothing

Change the optimization sense of the model `m` to `sense`. `sense` must be
`:min` or `:max`.
"""
function change_objective_sense! end

"""
    get_constant_objective(m)::Float64

Return the constant (i.e. offset) component of the objective.
"""
get_constant_objective(m::LinQuadOptimizer) = m.objective_constant

"""
    get_linear_objective!(m, x::Vector{Float64})

Get the linear coefficients of the objective and store
in `x`.
"""
function get_linear_objective! end

"""
    get_quadratic_terms_objective(m)::SparseMatrixCSC{Float64,Int64}

Get quadratic terms of the objective function returned in sparse CSC format.
"""
function get_quadratic_terms_objective end

"""
    get_objectivesense(m)::MOI.OptimizationSense

Get the optimization sense of the model `m`.
"""
function get_objectivesense end

"""
    solve_mip_problem!(m)::Nothing

Solve a mixed-integer model `m`.
"""
function solve_mip_problem! end

"""
    solve_quadratic_problem!(m)::Nothing

Solve a model `m` with quadratic components.
"""
function solve_quadratic_problem! end

"""
    solve_linear_problem!(m)::Nothing

Solve a linear program `m`.
"""
function solve_linear_problem! end

"""
    get_variable_primal_solution!(m, x::Vector{Float64})

Get the primal solution for the variables in the model `m`, and
store in `x`. `x`must have one element for each variable.
"""
function get_variable_primal_solution! end

"""
    get_linear_primal_solution!(m, x::Vector{Float64})

Given a set of linear constraints `l <= a'x <= b` in the model `m`, get the
constraint primal `a'x` for each constraint, and store in `x`.
`x` must have one element for each linear constraint.
"""
function get_linear_primal_solution! end

"""
    get_quadratic_primal_solution!(m, x::Vector{Float64})

Given a set of quadratic constraints `l <= a'x + x'Qx <= b` in the model `m`,
get the constraint primal `a'x + x'Qx` for each constraint, and store in `x`.
`x` must have one element for each quadratic constraint.
"""
function get_quadratic_primal_solution! end

"""
    get_variable_dual_solution!(m, x::Vector{Float64})

Get the dual solution (reduced-costs) for the variables in the model `m`, and
store in `x`. `x`must have one element for each variable.
"""
function get_variable_dual_solution! end

"""
    get_linear_dual_solution!(m, x::Vector{Float64})

Get the dual solution for the linear constraints in the model `m`, and
store in `x`. `x`must have one element for each linear constraint.
"""
function get_linear_dual_solution! end

"""
    get_quadratic_dual_solution!(m, x::Vector{Float64})

Get the dual solution for the quadratic constraints in the model `m`, and
store in `x`. `x`must have one element for each quadratic constraint.
"""
function get_quadratic_dual_solution! end

"""
    get_objective_value(m)

Get the objective value of the solved model `m`.
"""
function get_objective_value end

"""
    get_objective_bound(m)

Get the objective bound of the model `m`.
"""
function get_objective_bound end

"""
    get_relative_mip_gap(m)

Get the relative MIP gap of the solved model `m`.
"""
function get_relative_mip_gap end

"""
    get_iteration_count(m)

Get the number of simplex iterations performed during the most recent
optimization of the model `m`.
"""
function get_iteration_count end

"""
    get_barrier_iterations(m)

Get the number of barrier iterations performed during the most recent
optimization of the model `m`.
"""
function get_barrier_iterations end

"""
    get_node_count(m)

Get the number of branch-and-cut nodes expolored during the most recent
optimization of the model `m`.
"""
function get_node_count end

"""
    get_farkas_dual!(m, x::Vector{Float64})

Get the farkas dual (certificate of primal infeasiblility) for the linear
constraints in the model `m`, and store in `x`. `x`must have one element for
each linear constraint.
"""
function get_farkas_dual! end

"""
    get_farkas_dual_bounds!(model, dest::Vector{Float64})

Get the farkas dual (certificate of primal infeasibility) for the variable
bounds in the model `model`, and store in `dest`. `dest` must have one element
for each variable.

Since most solvers do not have this feature, this function has a default
fallback that does nothing.
"""
get_farkas_dual_bounds!(model::LinQuadOptimizer, dest) = nothing

"""
    get_unbounded_ray!(m, x::Vector{Float64})

Get the unbounded ray (certificate of dual infeasiblility) for the linear
constraints in the model `m`, and store in `x`. `x`must have one element for
each variable.
"""
function get_unbounded_ray! end

"""
    get_termination_status(m)

Get the termination status of the model `m`.
"""
function get_termination_status end

"""
    get_primal_status(m)

Get the primal status of the model `m`.
"""
function get_primal_status end

"""
    get_dual_status(m)

Get the dual status of the model `m`.
"""
function get_dual_status end

"""
    get_number_variables(m)::Int

Get the number of variables in the model `m`.
"""
function get_number_variables end

"""
    add_variables!(m, n::Int)::Nothing

Add `n` new variables to the model `m`.
"""
function add_variables! end

"""
    delete_variables!(m, start_col::Int, end_col::Int)::Nothing

Delete the columns `start_col`, `start_col+1`, ..., `end_col` from the model `m`.
"""
function delete_variables! end

"""
    add_mip_starts!(model::M, cols::Vector{Int}, x::Vector{Float64})::Nothing

Add a primal start `x` for the variables in the columns `cols` of `model`.

Note that if this method is implemented, solvers of type `M` must also declare
that they support VariablePrimalStarts by overloading the following method:

    function MOI.supports(model::M,
                          ::MOI.VariablePrimalStart,
                          ::Type{MOI.VariableIndex})
        return true
    end
"""
function add_mip_starts! end