text:
using Images, ImageDraw, ImageMagick
#Additional required packages to run this example.
#Pkg.add("Images")
#Pkg.add("ImageMagick")

using AprilTags

#-------------------------------------------------------------------------------
## Example AprilTags.jl detection from an image using the default setup

# Simple method to show the image with the tags
function showImage(image, tags, K)
    # Convert image to RGB
    imageCol = RGB.(image)
    #traw color box on tag corners
    foreach(tag->drawTagBox!(imageCol, tag, width = 2, drawReticle = false), tags)
    foreach(tag->drawTagAxes!(imageCol,tag, K), tags)
    imageCol
end

fx = 524.040
fy = 524.040
cy = 319.254
cx = 251.227
K = [-fx 0  cx;
      0 fy cy]

detector = nothing
try
    # Create default detector
    global detector = AprilTagDetector()

    # 1. Run against a file
    global image = load(dirname(Base.source_path()) *"/../data/tagtest.jpg");
    global tags = detector(image)
    showImage(image, tags, K)

    # 2. Run against an image from memory
    file = open(dirname(Base.source_path()) *"/../data/tagtest.jpg")
    imgBytes = read(file)
    close(file)
    # Here's where we pretend it came from a stream
    image = readblob(imgBytes) #ImageMagick function
    tags = detector(image)
    showImage(image, tags, K)
finally
    ## free memmory
    freeDetector!(detector)
end

##-------------------------------------------------------------------------------
detections = nothing
td = nothing
tf = nothing

try
    # 3. Use low-level methods and direct access to wrapper
    global image = load(dirname(Base.source_path()) *"/../data/tagtest.jpg")
    # Create april tag detector
    global td = apriltag_detector_create()
    # Create tag family
    global tf = tag36h11_create()
    # add family to detector
    apriltag_detector_add_family(td, tf)
    # create image8 object for april tags
    image8 = convert(AprilTags.image_u8, image)
    # run detector on image
    global detections =  apriltag_detector_detect(td, image8)
    # copy detections
    global tags = getTagDetections(detections)
    # Reading homography of tag 1 (deepcopy since memory is destoyed by c)
    voidpointertoH = Base.unsafe_convert(Ptr{Nothing}, tags[1].H)
    # pointer to H matrix
    nrows = unsafe_load(Ptr{UInt32}(voidpointertoH),1)
    ncols = unsafe_load(Ptr{UInt32}(voidpointertoH),2)
    H = deepcopy(unsafe_wrap(Array, Ptr{Cdouble}(voidpointertoH+8), (3,3)))

finally
    # Cleanup: free the detector and tag family when done.
    apriltag_detections_destroy(detections)
    apriltag_detector_destroy(td)
    tag36h11_destroy(tf)
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

@testset "Indicator by SOS1" begin
    # linear problem with indicator constraints
    # max  2x1 + 3x2
    # s.t. x1 + x2 <= 10
    #      z1 == 1 ==> x2 <= 8
    #      z2 == 1 ==> x1 + x2 == 9
    #      z3 == 1 ==> x1 >= 5
    #      z1 + z2 >= 1
    model = MOIU.MockOptimizer(MOIU.Model{Float64}());
    config = MOIT.TestConfig()

    x1 = MOI.add_variable(model)
    x2 = MOI.add_variable(model)
    z1 = MOI.add_variable(model)
    z2 = MOI.add_variable(model)
    z3 = MOI.add_variable(model)

    vc2 = MOI.add_constraint(model, z2, MOI.ZeroOne())
    vc1 = MOI.add_constraint(model, z1, MOI.ZeroOne())
    vc3 = MOI.add_constraint(model, z3, MOI.ZeroOne())
    f1 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0],
    )
    iset1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(8.0))

    f2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
        ],
        [0.0, 0.0]
    )
    iset2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(9.0))

    f3 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z3)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x1)),
        ],
        [0.0, 0.0]
    )
    iset3 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(5.0))

    BT1 = MOIB.Constraint.concrete_bridge_type(MOIB.Constraint.IndicatorSOS1Bridge{Float64}, typeof(f1), typeof(iset1))
    bridge1 = MOIB.Constraint.bridge_constraint(BT1, model, f1, iset1)
    @test BT1 <: MOIB.Constraint.IndicatorSOS1Bridge{Float64, <:MOI.LessThan, <:MOI.ConstraintIndex}
    @test bridge1 isa BT1

    BT2 = MOIB.Constraint.concrete_bridge_type(MOIB.Constraint.IndicatorSOS1Bridge{Float64}, typeof(f2), typeof(iset2))
    bridge2 = MOIB.Constraint.bridge_constraint(BT2, model, f2, iset2)
    @test BT2 <: MOIB.Constraint.IndicatorSOS1Bridge{Float64, <:MOI.EqualTo, Nothing}
    @test bridge2 isa BT2

    BT3 = MOIB.Constraint.concrete_bridge_type(MOIB.Constraint.IndicatorSOS1Bridge{Float64}, typeof(f3), typeof(iset3))
    bridge3 = MOIB.Constraint.bridge_constraint(BT3, model, f3, iset3)
    @test BT3 <: MOIB.Constraint.IndicatorSOS1Bridge{Float64, <:MOI.GreaterThan, <:MOI.ConstraintIndex}
    @test bridge3 isa BT3

    w1 = bridge1.w_variable_index
    @test MOI.get(model, MOI.ConstraintFunction(), bridge1.bound_constraint_index) == MOI.SingleVariable(w1)
    @test MOI.get(model, MOI.ConstraintSet(), bridge1.bound_constraint_index) == MOI.LessThan(0.0)
    @test MOI.get(model, MOI.ConstraintFunction(), bridge1.sos_constraint_index) == MOI.VectorOfVariables([w1, z1])
    lin_cons1 = bridge1.linear_constraint_index
    lin_func1 = MOI.get(model, MOI.ConstraintFunction(), lin_cons1)
    @test lin_func1 ≈ MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x2), MOI.ScalarAffineTerm(1.0, w1)], 0.0)
    @test MOI.get(model, MOI.ConstraintSet(), lin_cons1) == MOI.LessThan(8.0)

    w2 = bridge2.w_variable_index
    @test bridge2.bound_constraint_index === nothing
    @test MOI.get(model, MOI.ConstraintFunction(), bridge2.sos_constraint_index) == MOI.VectorOfVariables([w2, z2])
    lin_cons2 = bridge2.linear_constraint_index
    lin_func2 = MOI.get(model, MOI.ConstraintFunction(), lin_cons2)
    @test lin_func2 ≈ MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x1), MOI.ScalarAffineTerm(1.0, x2), MOI.ScalarAffineTerm(1.0, w2)], 0.0)
    @test MOI.get(model, MOI.ConstraintSet(), lin_cons2) == MOI.EqualTo(9.0)

    w3 = bridge3.w_variable_index
    @test MOI.get(model, MOI.ConstraintFunction(), bridge3.bound_constraint_index) == MOI.SingleVariable(w3)
    @test MOI.get(model, MOI.ConstraintSet(), bridge3.bound_constraint_index) == MOI.GreaterThan(0.0)
    @test MOI.get(model, MOI.ConstraintFunction(), bridge3.sos_constraint_index) == MOI.VectorOfVariables([w3, z3])
    lin_cons3 = bridge3.linear_constraint_index
    lin_func3 = MOI.get(model, MOI.ConstraintFunction(), lin_cons3)
    @test lin_func3 ≈ MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x1), MOI.ScalarAffineTerm(1.0, w3)], 0.0)
    @test MOI.get(model, MOI.ConstraintSet(), lin_cons3) == MOI.GreaterThan(5.0)

    ## MOI.get on bridge
    @test MOI.get(model, MOI.ConstraintSet(), bridge3) == iset3
    @test MOI.get(model, MOI.ConstraintFunction(), bridge3) ≈ f3
end

@testset "Basic constraint test" begin
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    config = MOIT.TestConfig()
    for BC in [MOI.LessThan{Float64}, MOI.GreaterThan{Float64}]
        bridged_mock = MOIB.Constraint.IndicatortoSOS1{Float64, BC, MOI.ConstraintIndex{MOI.SingleVariable, BC}}(mock)
        MOIT.basic_constraint_tests(bridged_mock, config,
            include=[
                (MOI.VectorAffineFunction{Float64}, MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, BC}),
            ],
        )
    end
end

@testset "Model equality" begin
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    config = MOIT.TestConfig()
    bridged_mock = MOIB.Constraint.IndicatortoSOS1{Float64, MOI.LessThan{Float64}, MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}}(mock)
    (z, bin_cons) = MOI.add_constrained_variable(bridged_mock, MOI.ZeroOne())
    MOI.set(bridged_mock, MOI.VariableName(), z, "z")
    MOI.set(mock, MOI.ConstraintName(), bin_cons, "bin_cons")
    x = MOI.add_variable(bridged_mock)
    MOI.set(bridged_mock, MOI.VariableName(), x, "x")
    var_names = ["z", "x", "w"]
    f = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x)),
        ],
        [0.0, 0.0]
    )
    iset = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(8.0))
    ci = MOI.add_constraint(bridged_mock, f, iset)
    @test length(MOI.get(mock, MOI.ListOfVariableIndices())) == 3
    MOI.set(mock, MOI.VariableName(), MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
    sos_cons_list = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorOfVariables, MOI.SOS1{Float64}}())
    @test length(sos_cons_list) == 1
    MOI.set(mock, MOI.ConstraintName(), sos_cons_list[1], "sos1")
    single_less_cons = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.LessThan{Float64}}())
    @test length(single_less_cons) == 1
    MOI.set(mock, MOI.ConstraintName(), single_less_cons[1], "wless")

    inequality_list = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
    @test length(inequality_list) == 1
    MOI.set(mock, MOI.ConstraintName(), inequality_list[1], "ineq")
    MOI.set(mock, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(z))
    MOI.set(mock, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    s = """
    variables: z, x, w
    sos1: [w, z] in MathOptInterface.SOS1([0.4, 0.6])
    wless: w <= 0.0
    ineq: x + w <= 8.0
    bin_cons: z in MathOptInterface.ZeroOne()
    maxobjective: z
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, s)
    MOIU.test_models_equal(mock, model, var_names, ["sos1", "wless", "ineq", "bin_cons"])

    test_delete_bridge(
        bridged_mock, ci, 2, (
            (MOI.SingleVariable, MOI.ZeroOne, 1),
            (MOI.SingleVariable, MathOptInterface.LessThan{Float64}, 0),
            (MOI.VectorOfVariables, MathOptInterface.SOS1{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
        ), used_bridges = 1, num_bridged = 1,
    )

    model = MOIU.Model{Float64}()
    sbridged = """
    variables: x, z
    bin_cons: z in MathOptInterface.ZeroOne()
    maxobjective: z
    """
    MOIU.loadfromstring!(model, sbridged)
    MOIU.test_models_equal(bridged_mock , model, ["z", "x"], ["bin_cons"])
end

@testset "Getting primal attributes" begin
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    config = MOIT.TestConfig()
    (z, bin_cons) = MOI.add_constrained_variable(mock, MOI.ZeroOne())
    x = MOI.add_variable(mock)
    f = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, z)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x)),
        ],
        [0.0, 0.0]
    )
    iset = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(8.0))
    BT = MOIB.Constraint.concrete_bridge_type(MOIB.Constraint.IndicatorSOS1Bridge{Float64}, typeof(f), typeof(iset))
    bridge1 = MOIB.Constraint.bridge_constraint(BT, mock, f, iset)
    # w value should be defaulted to 0
    MOI.set(mock, MOI.VariablePrimalStart(), bridge1.w_variable_index, 0.0)
    affine_value = 6.0
    MOI.set(mock, MOI.ConstraintPrimalStart(), bridge1, [1.0, affine_value])
    @test MOI.get(mock, MOI.VariablePrimalStart(), z) ≈ 1.0
    @test MOI.get(mock, MOI.ConstraintPrimalStart(), bridge1.linear_constraint_index) ≈ 6.0

    # after setting the w value
    w_value = 3.0
    MOI.set(mock, MOI.VariablePrimalStart(), bridge1.w_variable_index, w_value)
    # linear function should not move
    @test all(MOI.get(mock, MOI.ConstraintPrimalStart(), bridge1) .≈ (1.0, affine_value - w_value))
    iseteq = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(8.0))
    BT = MOIB.Constraint.concrete_bridge_type(MOIB.Constraint.IndicatorSOS1Bridge{Float64}, typeof(f), typeof(iseteq))
    bridge_eq = MOIB.Constraint.bridge_constraint(BT, mock, f, iseteq)
    @test MOI.get(bridge_eq, MOI.NumberOfConstraints{MOI.SingleVariable, MOI.EqualTo{Float64}}()) == 0
    @test MOI.get(bridge_eq, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}()) == 1
    @test isempty(MOI.get(bridge_eq, MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.EqualTo{Float64}}()))
    @test MOI.supports(mock, MOI.ConstraintPrimalStart(), MOIB.Constraint.IndicatorSOS1Bridge{Float64})
    @test MOI.supports(mock, MOI.ConstraintPrimalStart(), MOIB.Constraint.IndicatorSOS1Bridge)

    # VariablePrimal
    MOI.set(mock, MOI.VariablePrimal(), bridge1.w_variable_index, 33.0)
    z_value = 1.0
    MOI.set(mock, MOI.VariablePrimal(), bridge1.z_variable_index, z_value)
    MOI.set(mock, MOI.VariablePrimal(), x, affine_value)

    # linear function should not move
    @test all(MOI.get(mock, MOI.ConstraintPrimal(), bridge1) .≈ (1.0, affine_value))
end

ext: jl


text:
macro umode_labesl()
	:(Dict(:sy=>"sync",:fd=>"delay",:md=>"mean_delay",:sd=>"strob",:smd=>"strob_delay",:ga=>"general_async",:roa=>"random_order_async",:ssy=>"semi_sync"))
end

macro update_schemes()
	:(Dict(:sy=>:sy,:fd=>fixed_delay!,:md=>smoothed_fixed_delay!,:st=>strob!,:mst=>smoothed_strob!,:ga=>:ga,:roa=>:roa,:ssy=>semi_sync!))
end

function init_net!(net::Net2,node_indices::Vector{Int}=collect(vertices(net.graph));update_mode::Symbol=:sy,t0_mode::Symbol=:custom,init_h::Symbol=:last,init_array::Vector{Int}=Int[],init_constraints...)
	allunique(node_indices) || error("Duplicate node indices")
	intersect(vertices(net.graph),node_indices)==node_indices || error("Invalid node indices")
	umodes = @update_schemes
	haskey(umodes,update_mode) || error("Unexisting update scheme. Available schemes are :",keys(umodes))
	nodes_to_constraint = Int[]
	max_buffer_tlength = size(net.sim_buffer)[2]
	if !isempty(init_constraints)
		ids = [string(id) for (id,state) in init_constraints]
		allunique(ids) || error("Duplicate constraint specification")
		intersect(ids,(sigma.id for sigma in net.nodes))==ids || error("Trying to set unexisting nodes:",ids)
		all(i->i<:Union{Int,Bool,Vector{Int},Vector{Bool}},(typeof(state) for (id,state) in init_constraints)) || error("Wrong format to specify state")
		push!(nodes_to_constraint,(findfirst([sigma.id for sigma in net.nodes],id) for id in ids)...)
		all(map((x,y)->typeof(x) <: Vector ? intersect(x,y)==x : intersect([x],y)==[x],(state for (id,state) in init_constraints),(sigma.state_range for sigma in net.nodes[nodes_to_constraint]))) || error("Trying to set node to a forbidden state")
		node_indices = setdiff(node_indices,nodes_to_constraint)
		states = Vector{Vector{Int}}(length(init_constraints))
		for (ni,nsigma) in enumerate(nodes_to_constraint)
			state = init_constraints[ni][2]
			states[ni] = isempty(state) ? [net.nodes[nsigma].def_state] : (typeof(state)<:Vector ? state : [state])
		end
	end
	if max_buffer_tlength > 1
		if init_h==:default
			for ni in vcat(node_indices,nodes_to_constraint)
				net.sim_buffer[ni,1:end-1] = net.nodes[ni].def_state
			end
		elseif init_h==:rand
			for ni in node_indices
				net.sim_buffer[ni,1:end-1] = rand(net.nodes[ni].state_range,max_buffer_tlength-1)
			end
		elseif init_h!=:last
			error("Invalid initialization flag. Available options are\n:default for default, :rand for random")
		end
	end
	if t0_mode==:default
		net.sim_buffer[node_indices,end] = [sigma.def_state for sigma in net.nodes[node_indices]]
	elseif t0_mode==:custom
		if !isempty(init_array)
			length(init_array)==length(node_indices)||error("Initial configuration vector must have the same length of nodes vector")
			for (n,sigma) in enumerate(net.nodes[node_indices])
				if init_array[n] in sigma.state_range
					net.sim_buffer[node_indices[n],end] = init_array[n]
				else
					error("Trying to set node ",sigma.id," to a forbidden state, value ",init_array[n]," is out of possible states range ",sigma.state_range)
				end
			end
		end
	else
		error("Invalid initialization flag. Available options are\n:default for default, :rand for random :custom start from last available configuration")
	end
	for (ni,nsigma) in enumerate(nodes_to_constraint)
		state = states[ni]
		if length(state)<=max_buffer_tlength
			net.sim_buffer[ni,end-length(state)+1:end] = collect(state)
		else
			net.sim_buffer[ni,:] = collect(state[end-max_buffer_tlength+1:end])
		end
	end

	in_state_vals = nonzeros(net.in_states)
	# mu_in_state_vals = nonzeros(net.mu_in_states)
	nsigma_ins = rowvals(net.in_states)
	# tau_vals = nonzeros(net.in_taus)
	# if update_mode == :sy
	# 	for i in vertices(net.graph)
	# 		# j_in = nzrange(net.in_states,i)
	# 		# mu_in_state_vals[j_in] = in_state_vals[j_in] = net.sim_buffer[nsigma_ins[j_in],end]
	# 	end
	if update_mode == :fd || update_mode == :st
		for i in vertices(net.graph)
			for j in nzrange(net.in_states,i)
				in_state_vals[j] = net.sim_buffer[nsigma_ins[j],end-net.in_taus[nsigma_ins[j],i]]
			end
		end
	elseif update_mode == :md || update_mode == :mst
		mu_in_state_vals = nonzeros(net.mu_in_states)
		for i in vertices(net.graph)
			for j in nzrange(net.in_states,i)
				nsigma = nsigma_ins[j]
				tau = net.in_taus[nsigma,i]
				sigma_in = net.nodes[nsigma]
				if tau==0
					mu_in_state_vals[j] = in_state_vals[j] = net.sim_buffer[nsigma,end]
				else
					mu_in_state_vals[j] = mean_sigmak = mean(net.sim_buffer[nsigma,end-tau:end])
					in_state_vals[j] = sigma_in.state_range[findlast(mean_sigmak .>= sigma_in.thresholds)+1]
				end
			end
		end
	end
	return net.sim_buffer[:,end]
end

function rand_init_net!(net::Net2,node_indices::Vector{Int}=collect(vertices(net.graph));update_mode::Symbol=:sy,init_h::Symbol=:rand,init_constraints...)
	allunique(node_indices) || error("Duplicate node indices")
	intersect(vertices(net.graph),node_indices)==node_indices || error("Invalid node indices")
	umodes = @update_schemes
	haskey(umodes,update_mode) || error("Unexisting update scheme. Available schemes are :",keys(umodes))
	nodes_to_constraint = Int[]
	max_buffer_tlength = size(net.sim_buffer)[2]
	if !isempty(init_constraints)
		ids = [string(id) for (id,state) in init_constraints]
		allunique(ids) || error("Duplicate constraint specification")
		intersect(ids,(sigma.id for sigma in net.nodes))==ids || error("Trying to set unexisting node $(setdiff(ids,(sigma.id for sigma in net.nodes)))")
		all(i->i<:Union{Int,Bool,Vector{Int},Vector{Bool}},(typeof(state) for (id,state) in init_constraints)) || error("Wrong format to specify state")
		push!(nodes_to_constraint,(findfirst([sigma.id for sigma in net.nodes],id) for id in ids)...)
		all(map((x,y)->typeof(x) <: Vector ? intersect(x,y)==x : intersect([x],y)==[x],(state for (id,state) in init_constraints),(sigma.state_range for sigma in net.nodes[nodes_to_constraint]))) || error("Trying to set node to a forbidden state")
		node_indices = setdiff(node_indices,nodes_to_constraint)
		states = Vector{Vector{Int}}(length(init_constraints))
		for (ni,nsigma) in enumerate(nodes_to_constraint)
			state = init_constraints[ni][2]
			states[ni] = isempty(state) ? [net.nodes[nsigma].def_state] : unique(state)
			net.sim_buffer[nsigma,end] = length(state)>1 ? rand(state) : state[1]
		end
	end
	if max_buffer_tlength > 1
		if init_h==:default
			for ni in vcat(node_indices,nodes_to_constraint)
				net.sim_buffer[ni,1:end-1] = net.nodes[ni].def_state
			end
		elseif init_h==:rand
			for ni in node_indices
				net.sim_buffer[ni,1:end-1] = rand(net.nodes[ni].state_range,max_buffer_tlength-1)
			end
			for (n,ni) in enumerate(nodes_to_constraint)
				net.sim_buffer[ni,1:end-1] = rand(states[n],max_buffer_tlength-1)
			end
		elseif init_h!=:last
			error("Invalid initialization flag. Available options are\n:default for default, :rand for random")
		end
	end
	net.sim_buffer[node_indices,end] = [rand(sigma.state_range) for sigma in net.nodes[node_indices]]
	in_state_vals = nonzeros(net.in_states)
	nsigma_ins = rowvals(net.in_states)
	# tau_vals = nonzeros(net.in_taus)
	if update_mode == :fd || update_mode == :st
		for i in vertices(net.graph)
			for j in nzrange(net.in_states,i)
				in_state_vals[j] = net.sim_buffer[nsigma_ins[j],end-net.in_taus[nsigma_ins[j],i]]
			end
		end
	elseif update_mode == :md || update_mode == :mst
		mu_in_state_vals = nonzeros(net.mu_in_states)
		for i in vertices(net.graph)
			for j in nzrange(net.in_states,i)
				nsigma = nsigma_ins[j]
				tau = net.in_taus[nsigma,i]
				sigma_in = net.nodes[nsigma]
				if tau==0
					mu_in_state_vals[j] = in_state_vals[j] = net.sim_buffer[nsigma,end]
				else
					mu_in_state_vals[j] = mean_sigmak = mean(net.sim_buffer[nsigma,end-tau:end])
					in_state_vals[j] = sigma_in.state_range[findlast(mean_sigmak .>= sigma_in.thresholds)+1]
				end
			end
		end
	end
	# return net.sim_buffer[:,end]
	return nothing
end

macro updating_scheme(x)
	update_ex0 = quote
	end
	if eval(x) == :sy
		update_ex1 = quote
			setdiff(sigma_indices,forced_indices)
		end
		update_ex2 = quote
			sim[i,t+1] = net.nodes[i].rule[sim[nsigma_ins[nzrange(net.in_states,i)],t]]
		end
	elseif eval(x) == :ga
		update_ex1 = quote
			rand(sigma_indices)
		end
		update_ex2 = quote
			unchanged = setdiff(sigma_indices,forced_indices)
			if !(i in forced_indices)
				sim[i,t+1] = net.nodes[i].rule[sim[nsigma_ins[nzrange(net.in_states,i)],t]]
				unchanged = setdiff(unchanged,i)
			end
			# @show unchanged
			sim[unchanged,t+1] = sim[unchanged,t]
		end
	elseif eval(x) ==:roa
		update_ex0 = quote
			sim[:,t+1] = sim[:,t]
		end
		update_ex1 = quote
			shuffle(sigma_indices)
		end
		update_ex2 = quote
			if !(i in forced_indices)
				sim[i,t+1] = net.nodes[i].rule[sim[nsigma_ins[nzrange(net.in_states,i)],t+1]]
			end
		end
	end
	return esc(quote
		nsigma_ins = rowvals(net.in_states)
		sigma_noise_indices = [ni for (ni,eta) in sigma_noises]
		sigma_indices = 1:length(net.nodes)
		for t in i_t0:i_t0+steps-1
			$update_ex0
			forced_indices = Int[]
			# Updating forced nodes
			for (nsigma,rhythm_values) in nsigma_ext_forcing
				rhythm,forced_values = rhythm_values
				index_rhythm = mod(t_init+t-i_t0,length(rhythm))+1
				if rhythm[index_rhythm]
					#println("At t=",t_init+t+1-i_t0,",trying to fix node index ",nsigma," on rhythm index ",index_rhythm," to forced_value ",forced_values[1])
					sim[nsigma,t+1] = forced_values[1]
					push!(forced_indices,nsigma)
					nsigma_ext_forcing[nsigma][2] = circshift(forced_values,-1)
				end
			end
			# Apply the original rule except for the forced nodes
			for i in $update_ex1
				$update_ex2
				# if i in forced_indices
				# 	println("At t=",t_init+t+1-i_t0,",node index ",i," changed from ",sim[i,t]," to ",sim[i,t+1])
				# end
				# Add noise where specified, expect for those forced nodes
				if i in sigma_noise_indices
					eta = sigma_noises[findfirst(sigma_noise_indices,i)][2]
					true_value = sim[i,t+1]
					sim[i,t+1] = rand()<=eta ? rand(filter(x -> x!=true_value,net.nodes[i].state_range)) : true_value
				end
				# sim[i,t+1] = net.nodes[i].rule[sim[nsigma_ins[nzrange(net.in_states,i)],t]]
			end
		end
	end)
end

function evolve_net!(net::Net2{T},steps::Int;keep::Bool=true,update_mode::Symbol=:sy,sim::Array{T,2}=Array{T}(nv(net.graph),steps + 1 + (update_mode == :sy ? 0 : maximum(net.in_taus))),t_init::Int=0,noise_vector::Vector{Pair{Symbol,Float64}}=Pair{Symbol,Float64}[],forcing_rhythms::Dict{Symbol,Tuple{Vector{Bool},Vector{Int}}}=Dict{Symbol,Tuple{Vector{Bool},Vector{Int}}}()) where {T<:state_type}
	umodes = @update_schemes
	haskey(umodes,update_mode) || error("Unexisting update scheme. Available schemes are :",keys(umodes))
	steps>=0 || error("Steps must be a positive integer")
	hist_size = 1 + (update_mode == :sy ? 0 : maximum(net.in_taus))
	# Initialize simulation table (from t_init - tau0 to t_init)
	size(sim) == (nv(net.graph),steps + hist_size) || error("Array for storing results has incorrect dimensions ",size(sim), ". Expected array size is ",(nv(net.graph),steps + hist_size))
	if update_mode == :sy
		sim[:,1] = net.sim_buffer[:,end]
	else
		sim[:,1:hist_size] = net.sim_buffer[:]
	end
	sigma_noises = collect(findfirst(x->x.id==string(sigma),net.nodes)=>eta for (sigma,eta) in noise_vector)
	nsigma_ext_forcing = Dict(findfirst(x->x.id==string(sigma_symbol),net.nodes)=>collect(rhythm_states) for (sigma_symbol,rhythm_states) in forcing_rhythms)
	i_t0 = hist_size
	if update_mode == :sy
		@updating_scheme :sy
	elseif update_mode == :ga
		for (nsigma,rhythm_states) in nsigma_ext_forcing
			nsigma_ext_forcing[nsigma][1] = repeat(nsigma_ext_forcing[nsigma][1],inner=[length(net.nodes)])
			nsigma_ext_forcing[nsigma][2] = repeat(nsigma_ext_forcing[nsigma][2],inner=[length(net.nodes)])
		end
		@show nsigma_ext_forcing
		@updating_scheme :ga
	elseif update_mode == :roa
		@updating_scheme :roa
	else
		umodes[update_mode](net,sim,steps,hist_size,sigma_noises,sigma_ext_forcing,t_init)
	end
	buffer_size = size(net.sim_buffer)[2]
	if buffer_size == hist_size || steps > buffer_size
		net.sim_buffer[:] = sim[:,end-buffer_size+1:end]
	elseif steps <= hist_size
		net.sim_buffer[:,1:end-steps-hist_size] = net.sim_buffer[:,1+steps+hist_size:end]
		net.sim_buffer[:,end-steps-hist_size+1:end] = sim[:]
	else
		@show steps, hist_size,buffer_size
		error("Wrong dimensions")
	end
	if keep
		return sim
	else
		return sim[:,end-hist_size+1:end]
	end
end

function semi_sync!(net::Net2,sim::sim_type,steps::Int,i_t0::Int,x...)
	nsigma_ins = rowvals(net.in_states)
	Nsize = length(net.nodes)
	buff = deepcopy(sim[:,i_t0])
	for t in i_t0:i_t0+steps-1
		n = 0
		for s in 1:Nsize
			s_indices = find(x->x.sync_order==s,net.nodes)
			for i in s_indices
				sim[i,t+1] = net.nodes[i].rule[buff[nsigma_ins[nzrange(net.in_states,i)]]]
				# sim[i,t+1] = net.nodes[i].rule[sim[nsigma_ins[nzrange(net.in_states,i)],t]]
				n += 1
			end
			buff[s_indices] = sim[s_indices,t+1]
			if n==Nsize
				break
			end
			# sim[:,t+1] = buff[:]
		end
	end
	return nothing
end

function fixed_delay!(net::Net2,sim::sim_type,steps::Int,i_t0::Int,x...)
	nsigma_ins = rowvals(net.in_states)
	# tau_ins = nonzeros(net.in_taus)
	# in_state_vals = nonzeros(net.in_states)
	# mu_in_state_vals = nonzeros(net.mu_in_states)
	for t in i_t0:i_t0+steps-1
		for (i,sigma) in enumerate(net.nodes)
			j_in = nzrange(net.in_states,i)
			input = [sim[nsigma_ins[nj],t-net.in_taus[nsigma_ins[nj],i]] for nj in j_in]
			sim[i,t+1] = sigma.rule[input]
			# mu_in_state_vals[j_in] = in_state_vals[j_in] = input
			# in_state_vals[j_in] = input
		end
	end
	return nothing
end

function smoothed_fixed_delay!(net::Net2,sim::sim_type,steps::Int,i_t0::Int,x...)
	nsigma_ins = rowvals(net.in_states)
	# tau_ins = nonzeros(net.in_taus)
	in_state_vals = nonzeros(net.in_states)
	mu_in_state_vals = nonzeros(net.mu_in_states)
	for t in i_t0:i_t0+steps-1
		for (i,sigma) in enumerate(net.nodes)
			j_in = nzrange(net.in_states,i)
			sim[i,t+1] = sigma.rule[in_state_vals[j_in]]
		end
		for (i,sigma) in enumerate(net.nodes)
			for j in nzrange(net.in_states,i)
				nsigma = nsigma_ins[j]
				tau = net.in_taus[nsigma,i]
				sigma_in = net.nodes[nsigma]
				if tau==0
					mu_in_state_vals[j] = in_state_vals[j] = sim[nsigma,t+1]
				else
					mu_in_state_vals[j] += (sim[nsigma,t+1]-sim[nsigma,t-tau])/(tau+1)
					in_state_vals[j] = sigma_in.state_range[findlast(mu_in_state_vals[j] .>= sigma_in.thresholds)+1]
				end
			end
		end
	end
	return nothing
end

function strob!(net::Net2,sim::sim_type,steps::Int,i_t0::Int,t_init::Int=0,x...)
	nsigma_ins = rowvals(net.in_states)
	tau_ins = nonzeros(net.in_taus)
	in_state_vals = nonzeros(net.in_states)
	for t in i_t0:i_t0+steps-1
		for (i,sigma) in enumerate(net.nodes)
			j_in = nzrange(net.in_states,i)
			for j in j_in
				if mod(t_init+t-i_t0+1,net.in_taus[nsigma_ins[j],i]+1)==0
					# sigma.sigma_in_state[nsigmak] = sim[sigmak,t-sigma.tau_in[nsigmak]]
					in_state_vals[j] = sim[nsigma_ins[j],t]
				end
			end
			sim[i,t+1] = sigma.rule[in_state_vals[j_in]]
		end
	end
	return nothing
end

function smoothed_strob!(net::Net2,sim::sim_type,steps::Int,i_t0::Int,t_init::Int=0,x...)
	nsigma_ins = rowvals(net.in_states)
	# tau_ins = nonzeros(net.in_taus)
	in_state_vals = nonzeros(net.in_states)
	mu_in_state_vals = nonzeros(net.mu_in_states)
	for t in i_t0:i_t0+steps-1
		for (i,sigma) in enumerate(net.nodes)
			j_in = nzrange(net.in_states,i)
			for j in j_in
				nsigma = nsigma_ins[j]
				tau = net.in_taus[nsigma,i]
				sigma_in = net.nodes[nsigma]
				if tau==0
					mu_in_state_vals[j] = in_state_vals[j] = sim[nsigma,t]
				elseif mod(t_init+t-i_t0+1,tau+1)==0
					mean_sigmak = mu_in_state_vals[j] = mean(sim[nsigma,t-tau:t])
					in_state_vals[j] = sigma_in.state_range[findlast(mean_sigmak.>= sigma_in.thresholds)+1]
				end
			end
			sim[i,t+1] = sigma.rule[in_state_vals[j_in]]
		end
	end
	return nothing
end

function obs_rand_cond!(net::Net2{T},steps::Int,obs_node::Union{AbstractString,Int,Vector{String},Vector{Int}},nrand_init::Int=prod(length(sigma.state_range) for sigma in net.nodes);base_dir::AbstractString="",variant_name::AbstractString="",file_tag::AbstractString="",update_mode::Symbol=:sy,noise_vector::Vector{Pair{Symbol,Float64}}=Pair{Symbol,Float64}[],init_h::Symbol=:default,forcing_rhythms::Dict{Symbol,Tuple{Vector{Bool},Vector{Int}}}=Dict{Symbol,Tuple{Vector{Bool},Vector{Int}}}(),init_constraints...) where {T<:state_type}
	Nsize = length(net.nodes)
	dump_to_file = false
	moment = now(Dates.UTC)
	time_string = string(Dates.value(moment))
	base_dir_tag = join([variant_name,file_tag,time_string],"_")
	results_dir = joinpath(base_dir,base_dir_tag)
	if base_dir != ""
		if isdir(base_dir)
			dump_to_file = true
			if isdir(results_dir)
				println("Resuls will be overwritten onto ")
			else
				mkpath(results_dir)
				println("Resuls will be saved into ")
			end
			println(results_dir," folder")
		else
			error(base_dir, " is not a valid directory to save s")
		end
	end
	index_node_obs = Any[]
	if typeof(obs_node) <: Vector
		length(obs_node)!=0 || error("You must enter at least 1 node id or index")
		index_node_obs = typeof(obs_node) <: Vector{Int} ? obs_node : collect(findfirst(sigma->sigma.id==x,net.nodes) for x in obs_node)
		for (i,x) in enumerate(obs_node)
			0 < index_node_obs[i] <= Nsize || error("Unexisting node: $x")
		end
	else
		index_node_obs = typeof(obs_node) <: Int ? obs_node : findfirst(sigma->sigma.id==obs_node,net.nodes)
		0 < index_node_obs <= Nsize || error("Unexisting node: $obs_node")
	end
	length_obs = length(index_node_obs)
	obs = length_obs >=2 ? Vector{Tuple{Vararg{T,length_obs}}}() : Vector{T}()
	state_space_size = prod(length(sigma.state_range) for sigma in net.nodes)
	if(state_space_size==0) && all(x->isinteger(x)&&x>0,length(sigma.state_range) for sigma in net.nodes)
		state_space_size = 10000000000
	end
	if nrand_init>state_space_size
        nrand_init = state_space_size
        println("Oversampling error: Random initial conditions bigger than possible number of configurations")
    end
	println("Calculating ",length_obs >=2 ? join((sigma.id for sigma in net.nodes[index_node_obs]),',') : net.nodes[index_node_obs].id ," dynamics over $nrand_init random initial conditions")
	r = 0
	hist_size = 1 + (update_mode == :sy ? 0 : maximum(net.in_taus))
	t_size = steps+hist_size
	sim_result = Array{T,2}(Nsize,t_size)
	sigma_ext_forcing = collect(sigma_symbol=>rhythm_values[2] for (sigma_symbol,rhythm_values) in forcing_rhythms)
	# init_constraints = vcat(init_constraints,sigma_ext_forcing)
	if dump_to_file
		#f = Vector{GZip.GZipStream}(length_obs)
		filenames = String[]
		raw_dir = joinpath(results_dir,"raw")
		isdir(raw_dir) || mkpath(raw_dir)
		for sigma_i in index_node_obs
			# push!(filenames,joinpath(results_dir,join([join([net.nodes[sigma_i].id,"T$t_size","r$nrand_init"],"_"),"bin","gz"],".")))
			push!(filenames,joinpath(raw_dir,join([net.nodes[sigma_i].id,"bin","gz"],".")))
			# f[i] = GZip.open(filename,"w")
			f = GZip.open(filenames[end],"w")
			close(f)
		end
		println("Succesfully created files for storing simulations")
		# for (i,sigma_i) in enumerate(index_node_obs)
		# 	filename = joinpath(results_dir,join([join([join([net.nodes[sigma_i].id,"T$t_size","r$nrand_init"],"_"),file_tag],""),"dat","gz"],"."))
		# 	f[i] = GZip.open(filename,"a")
		# end
		if nrand_init<state_space_size
			while r < nrand_init
				rand_init_net!(net;update_mode=update_mode,init_h=init_h,init_constraints...)
				evolve_net!(net,steps;keep=false,update_mode=update_mode,sim=sim_result,noise_vector=noise_vector,forcing_rhythms=forcing_rhythms)
				# nbytes = Int[]
				# for (fi,obs_i) in enumerate(index_node_obs)
				# 	#push!(nbytes,write(f[fi],map(Int8,sim_result[obs_i,:])))
				# 	writedlm(f[fi],map(Int8,sim_result[obs_i,:])',',')
				# end
				for (fi,sigma_i) in enumerate(index_node_obs)
					GZip.open(filenames[fi],"a") do filex
						write(filex,sim_result[sigma_i,:])
					end
					# close(f)
				end
				#println(nbytes," written")
				r += 1
				if nrand_init >= 10 && mod(r,round(Int,nrand_init/10))==0
					@printf("%0.2f%% progress, %d out of %d initial conditions explored\n", 100*r/nrand_init, r, nrand_init)
				end
			end
		else
			println("Exhaustive exploration of state space")
			node_indices = collect(vertices(net.graph))
			nodes_to_constraint = Int[]
			if !isempty(init_constraints)
				ids = [string(id) for (id,state) in init_constraints]
				allunique(ids) || error("Duplicate constraint specification")
				intersect(ids,(sigma.id for sigma in net.nodes))==ids || error("Trying to set unexisting node")
				all(i->i<:Union{Int,Bool,Vector{Int},Vector{Bool}},(typeof(state) for (id,state) in init_constraints)) || error("Wrong format to specify state")
				push!(nodes_to_constraint,(findfirst([sigma.id for sigma in net.nodes],id) for id in ids)...)
				all(map((x,y)->typeof(x) <: Vector ? intersect(x,y)==x : intersect([x],y)==[x],(state for (id,state) in init_constraints),(sigma.state_range for sigma in net.nodes[nodes_to_constraint]))) || error("Trying to set node to a forbidden state")
				node_indices = setdiff(node_indices,nodes_to_constraint)
				states = Vector{Vector{Int}}(length(init_constraints))
				for (ni,nsigma) in enumerate(nodes_to_constraint)
					state = init_constraints[ni][2]
					states[ni] = isempty(state) ? [net.nodes[nsigma].def_state] : unique(state)
				end
			end
			ranges = Vector{Vector{Int}}(nv(net.graph))
			for nsigma in node_indices
				ranges[nsigma] = net.nodes[nsigma].state_range
			end
			for (ni,nsigma) in enumerate(nodes_to_constraint)
				ranges[nsigma] = states[ni]
			end
			itstspace = product(ranges...)
			for leave in itstspace
				init_net!(net;update_mode=update_mode,init_h=init_h,init_array=collect(leave))
				evolve_net!(net,steps;keep=false,update_mode=update_mode,sim_result=sim_result,noise_vector=noise_vector,forcing_rhythms=forcing_rhythms)
				for (fi,sigma_i) in enumerate(index_node_obs)
					GZip.open(filenames[fi],"a") do filex
						write(filex,sim_result[sigma_i,:])
					end
					# close(f)
				end
				r += 1
				if nrand_init>=10 && mod(r,round(Int,nrand_init/10))==0
					@printf("%0.2f%% progress, %d out of %d initial conditions explored\n", 100*r/nrand_init, r, nrand_init)
				end
			end
			nrand_init = r
		end
		# status_files = zeros(Int,length_obs)
		# for (fi,sigma_i) in enumerate(index_node_obs)
		# 	status_files[fi] = close(f[fi])
		# end
		# if isempty(find(status_files))
		# 	println("Simulations were succesfully stored for nodes ", join(obs_node,", "," and "))
		# else
		# 	println("Simulations couldn't be stored for nodes ", join(obs_node[find(status_files)],", "," nor "))
		# end
		# return status_files
		println("\nWriting metadata files")
		open(joinpath(results_dir,join([base_dir_tag,"meta"],".")), "w") do f
			println(f, "Date = ",string(moment)) # Univeral Time
			println(f, "umode = ",@umode_labesl()[update_mode]) # tag for updating scheme used in the simulation
        	println(f, "t0 = ",1-hist_size) # What time corresponds to first time instant, e.g. 1st column of array element corresponds to t=0 in synchronous update scheme
			println(f, "T = ",steps) # Amount of iterations
			println(f, "tsize = ",hist_size + steps) # Whole time series length
			println(f, "r = ",nrand_init) # Total amount of realizations
			println(f, "\nSampling space of initial conditions") # For each variable, it's sampling space is printed out
			init_constraints_lab = [string(nsymbol) for (nsymbol,nvals) in init_constraints]
			init_constraints_states = [nvals for (nsymbol,nvals) in init_constraints]
			for node in net.nodes
				print(f,node.id,"\t:\t")
				if node.id in init_constraints_lab
					println(f,join(init_constraints_states[findfirst(x->x==node.id,init_constraints_lab)],","))
				else
					println(f,join(node.state_range,","))
				end
			end
		end
		return nothing
	else
		if nrand_init<state_space_size
			while r < nrand_init
				rand_init_net!(net;update_mode=update_mode,init_h=init_h,init_constraints...)
				evolve_net!(net,steps;keep=false,update_mode=update_mode,sim=sim_result,noise_vector=noise_vector,forcing_rhythms=forcing_rhythms)
				if length_obs==1
					push!(obs,sim_result[index_node_obs,:]...)
				else
					push!(obs,(tuple(sim_result[index_node_obs,ti]...) for ti in 1:t_size)...)
				end
				r += 1
				if nrand_init>=10 && mod(r,round(Int,nrand_init/10))==0
					@printf("%0.2f%% progress, %d out of %d initial conditions explored\n", 100*r/nrand_init, r, nrand_init)
				end
			end
		else
			println("Exhaustive exploration of state space")
			node_indices = collect(vertices(net.graph))
			nodes_to_constraint = Int[]
			if !isempty(init_constraints)
				ids = [string(id) for (id,state) in init_constraints]
				allunique(ids) || error("Duplicate constraint specification")
				intersect(ids,(sigma.id for sigma in net.nodes))==ids || error("Trying to set unexisting node")
				all(i->i<:Union{Int,Bool,Vector{Int},Vector{Bool}},(typeof(state) for (id,state) in init_constraints)) || error("Wrong format to specify state")
				push!(nodes_to_constraint,(findfirst([sigma.id for sigma in net.nodes],id) for id in ids)...)
				all(map((x,y)->typeof(x) <: Vector ? intersect(x,y)==x : intersect([x],y)==[x],(state for (id,state) in init_constraints),(sigma.state_range for sigma in net.nodes[nodes_to_constraint]))) || error("Trying to set node to a forbidden state")
				node_indices = setdiff(node_indices,nodes_to_constraint)
				states = Vector{Vector{Int}}(length(init_constraints))
				for (ni,nsigma) in enumerate(nodes_to_constraint)
					state = init_constraints[ni][2]
					states[ni] = isempty(state) ? [net.nodes[nsigma].def_state] : unique(state)
				end
			end
			ranges = Vector{Vector{Int}}(nv(net.graph))
			for nsigma in node_indices
				ranges[nsigma] = net.nodes[nsigma].state_range
			end
			for (ni,nsigma) in enumerate(nodes_to_constraint)
				ranges[nsigma] = states[ni]
			end
			itstspace = product(ranges...)
			for leave in itstspace
				init_net!(net;update_mode=update_mode,init_h=init_h,init_array=collect(leave))
				evolve_net!(net,steps;keep=false,update_mode=update_mode,sim_result=sim_result,noise_vector=noise_vector,forcing_rhythms=forcing_rhythms)
				if length_obs==1
					push!(obs,sim_result[index_node_obs,:]...)
				else
					push!(obs,(tuple(sim_result[index_node_obs,ti]...) for ti in 1:t_size)...)
				end
				#push!(obs,sim_result[index_node_obs,:]...)
				r += 1
				if nrand_init>=10 && mod(r,round(Int,nrand_init/10))==0
					@printf("%0.2f%% progress, %d out of %d initial conditions explored\n", 100*r/nrand_init, r, nrand_init)
				end
			end
			nrand_init = r
		end
		return reshape(obs,round(Int,length(obs)/nrand_init),nrand_init)
	end
end


text:
using LadeGCM

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

c = constants();
@testset "Configuration" begin
    @testset "Defaults" begin
        @test c.cₐ₀ ≈ 589.0
        @test c.cₜ₀ ≈ 1875.0
        @test c.cₘ₀ ≈ 900.0
        @test c.λ ≈ 1.8
        @test c.τ ≈ 4.0
        @test c.D ≈ 1.0
        @test c.r ≈ 12.5
        @test c.DT ≈ 0.0423
        @test c.B₀ ≈ 13.0
        @test c.BT ≈ 0.032
        @test c.w₀ ≈ 0.1
        @test c.wT ≈ 0.1
        @test c.QR ≈ 1.72
        @test c.NPP₀ ≈ 55.0
        @test c.KC ≈ 0.3
    end
    @testset "Overrides" begin
        co = constants(cₐ₀=2.0, cₜ₀=3.0, QR=4.0, DT=5.0);
        @test co.cₐ₀ ≈ 2.0 #Altered
        @test co.cₜ₀ ≈ 3.0 #Altered
        @test co.cₘ₀ ≈ 900.0
        @test co.λ ≈ 1.8
        @test co.τ ≈ 4.0
        @test co.D ≈ 1.0
        @test co.r ≈ 12.5
        @test co.DT ≈ 5.0 #Altered
        @test co.B₀ ≈ 13.0
        @test co.BT ≈ 0.032
        @test co.w₀ ≈ 0.1
        @test co.wT ≈ 0.1
        @test co.QR ≈ 4.0 #Altered
        @test co.NPP₀ ≈ 55.0
        @test co.KC ≈ 0.3
    end
end

@testset "Pathways" begin
    rcp26 = RCP3PD;
    rcp45 = RCP45;
    rcp6 = RCP6;
    rcp85 = RCP85;
    @testset "Types" begin
        @test typeof(rcp26) <: LadeGCM.Pathway
        @test isa(rcp26, LadeGCM.RCP3PDPathway)
        @test typeof(rcp45) <: LadeGCM.Pathway
        @test isa(rcp45, LadeGCM.RCP45Pathway)
        @test typeof(rcp6) <: LadeGCM.Pathway
        @test isa(rcp6, LadeGCM.RCP6Pathway)
        @test typeof(rcp85) <: LadeGCM.Pathway
        @test isa(rcp85, LadeGCM.RCP85Pathway)
    end
    @testset "Show" begin
        @test sprint(show, rcp26) == "+2.6 W/m² (Peak & Decline) Representative Concentration Pathway"
        @test sprint(show, rcp45) == "+4.5 W/m² Representative Concentration Pathway"
        @test sprint(show, rcp6) == "+6.0 W/m² Representative Concentration Pathway"
        @test sprint(show, rcp85) == "+8.5 W/m² Representative Concentration Pathway"
    end
    @testset "Summary" begin
        @test summary(rcp26) == "RCP3PD"
        @test summary(rcp45) == "RCP45"
        @test summary(rcp6) == "RCP6"
        @test summary(rcp85) == "RCP85"
    end
end

@testset "Emission Parameters" begin
    @testset "Loading data" begin
        data_frame = LadeGCM.load_pathway_data(RCP45);
        @test data_frame[:Year][1] == 1765
        @test data_frame[:FossilCO2][345] ≈ 3.19374
        @test data_frame[:OtherCO2][end] ≈ 0.0
    end
    @testset "Generate Continuous Data" begin
        (E, LUC) = LadeGCM.generate_emission_parameters(RCP45);
        @test E(1768.8) ≈ 0.003
        @test E(2100) ≈ 4.203
        @test LUC(1982.77) ≈ 1.217443742
        @test LUC(2019) ≈ 0.39778
    end
end

@testset "Solutions" begin
    # NOTE: These are going to have some numerical error in them,
    # so we round everything for stability.
    r6 = calculate(RCP6, saveat=1765:10:2100);
    r85 = calculate(RCP85, saveat=1765:10:2100);
@static if VERSION < v"0.7.0-DEV.2005"
    @test isapprox(round(r6.cₜ[13],0), 1863.0; atol = 2)
    @test isapprox(round(r6.cₘ[28],0), 930.0; atol = 2)
    @test isapprox(round(r6.cₐ[3],0), 590.0; atol = 2)
    @test isapprox(round(r6.ΔT[9],0), 0.0; atol = 2)
    @test isapprox(round(r6.Δcₜ[27],0), 2.0; atol = 2)
    @test isapprox(round(r6.Δcₘ[14],0), 0.0; atol = 2)
    @test isapprox(round(r6.ΔcM[23],0), 1.0; atol = 2)
    @test isapprox(round(r6.Δcₛ[16],0), 1.0; atol = 2)
    @test isapprox(round(r6.Δcₐ[31],0), 10.0; atol = 2)
    @test isapprox(round(r6.year[24],0), 1995.0; atol = 2)
    @test isapprox(round(r85.cₜ[30],0), 1991.0; atol = 2)
    @test isapprox(round(r85.cₘ[8],0), 900.0; atol = 2)
    @test isapprox(round(r85.cₐ[16],0), 614.0; atol = 2)
    @test isapprox(round(r85.ΔT[29],0), 2.0; atol = 2)
    @test isapprox(round(r85.Δcₜ[21],0), 0.0; atol = 2)
    @test isapprox(round(r85.Δcₘ[18],0), 0.0; atol = 2)
    @test isapprox(round(r85.ΔcM[4],0), 0.0; atol = 2)
    @test isapprox(round(r85.Δcₛ[26],0), 7.0; atol = 2)
    @test isapprox(round(r85.Δcₐ[10],0), 0.0; atol = 2)
    @test isapprox(round(r85.year[28],0), 2035.0; atol = 2)
else
    @test isapprox(round(r6.cₜ[13];digits=0), 1863.0; atol = 2)
    @test isapprox(round(r6.cₘ[28];digits=0), 930.0; atol = 2)
    @test isapprox(round(r6.cₐ[3];digits=0), 590.0; atol = 2)
    @test isapprox(round(r6.ΔT[9];digits=0), 0.0; atol = 2)
    @test isapprox(round(r6.Δcₜ[27];digits=0), 2.0; atol = 2)
    @test isapprox(round(r6.Δcₘ[14];digits=0), 0.0; atol = 2)
    @test isapprox(round(r6.ΔcM[23];digits=0), 1.0; atol = 2)
    @test isapprox(round(r6.Δcₛ[16];digits=0), 1.0; atol = 2)
    @test isapprox(round(r6.Δcₐ[31];digits=0), 10.0; atol = 2)
    @test isapprox(round(r6.year[24];digits=0), 1995.0; atol = 2)
    @test isapprox(round(r85.cₜ[30];digits=0), 1991.0; atol = 2)
    @test isapprox(round(r85.cₘ[8];digits=0), 900.0; atol = 2)
    @test isapprox(round(r85.cₐ[16];digits=0), 614.0; atol = 2)
    @test isapprox(round(r85.ΔT[29];digits=0), 2.0; atol = 2)
    @test isapprox(round(r85.Δcₜ[21];digits=0), 0.0; atol = 2)
    @test isapprox(round(r85.Δcₘ[18];digits=0), 0.0; atol = 2)
    @test isapprox(round(r85.ΔcM[4];digits=0), 0.0; atol = 2)
    @test isapprox(round(r85.Δcₛ[26];digits=0), 7.0; atol = 2)
    @test isapprox(round(r85.Δcₐ[10];digits=0), 0.0; atol = 2)
    @test isapprox(round(r85.year[28];digits=0), 2035.0; atol = 2)
end
end

max_stars_repo_name: Libbum/LadeGCM

max_issues_repo_name: Libbum/LadeGCM

max_forks_repo_name: Libbum/LadeGCM


text:
"""
    wolff!(model::Ising, avg_cluster_size::Int)

Implementation of the Wolff single cluster algorithm for the Ising model.

Arguments
* `model`           : Ising type
* `avg_cluster_size`: Average cluster size
"""
function wolff!(model::Ising, avg_cluster_size::Int)
    P = 1 - exp(-2 * model.params.β)
    n_clusters = Int(ceil(model.params.L^2 / avg_cluster_size))
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)
    indices = 1:length(model.σ)

    for t = 1:model.params.n_sweeps
        for n = 1:n_clusters
            idx = CartesianIndices(model.σ)[rand(model.rng, indices)]
            i, j = idx[1], idx[2]
            old_spin = model.σ[i, j]
            new_spin = -old_spin
            model.σ[i, j] = new_spin

            cluster = [(i, j)]
            for (i, j) in cluster
                for (k, l) in get_neighbor_indices(model, (i, j))
                    if model.σ[k, l] == old_spin && rand(model.rng) < P
                        model.σ[k, l] = new_spin
                        push!(cluster, (k, l))
                    end
                end
            end

            model.observables.avg_cluster_size += length(cluster)
        end

        t > t₀ && update_observables!(model)
    end

    model.observables.avg_cluster_size /= (model.params.n_sweeps * n_clusters)
    compute_observables_statistics!(model)
end

"""
    wolff!(model::Potts, avg_cluster_size::Int)

Implementation of the Wolff single cluster algorithm for the Potts model.
Arguments
* `model`           : Potts type
* `avg_cluster_size`: Average cluster size
"""
function wolff!(model::Potts, avg_cluster_size::Int)
    P = 1 - exp(-model.params.β)
    n_clusters = Int(ceil(model.params.L^2 / avg_cluster_size))
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)
    indices = 1:length(model.σ)

    for t = 1:model.params.n_sweeps
        for n = 1:n_clusters
            idx = CartesianIndices(model.σ)[rand(model.rng, indices)]
            i, j = idx[1], idx[2]
            old_spin = model.σ[i, j]
            new_spin = rand(model.rng, UnitRange{Int8}(1, model.params.q))
            model.σ[i, j] = new_spin

            cluster = [(i, j)]
            for (i, j) in cluster
                for (k, l) in get_neighbor_indices(model, (i, j))
                    if model.σ[k, l] == old_spin && rand(model.rng) < P
                        model.σ[k, l] = new_spin
                        push!(cluster, (k, l))
                    end
                end
            end

            model.observables.avg_cluster_size += length(cluster)
        end

        t > t₀ && update_observables!(model)
    end

    model.observables.avg_cluster_size /= (model.params.n_sweeps * n_clusters)
    compute_observables_statistics!(model)
end

"""
    wolff!(model::XY, avg_cluster_size::Int)

Implementation of the Wolff single cluster algorithm for the XY model.

Arguments
* `model`           : XY type
* `avg_cluster_size`: Average cluster size
"""
function wolff!(model::XY, avg_cluster_size::Int)
    P = 1 - exp(-2 * model.params.β)
    n_clusters = Int(ceil(model.params.L^2 / avg_cluster_size))
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)
    indices = 1:length(model.σ)

    for t = 1:model.params.n_sweeps
        for n = 1:n_clusters
            r = random_XYVector(model.rng)
            idx = CartesianIndices(model.σ)[rand(model.rng, indices)]
            i, j = idx[1], idx[2]
            flip_spin!(model, (i, j), r)

            stack = [(i, j)]
            cluster = Set(stack)
            while length(stack) > 0
                i, j = pop!(stack)
                spin = model.σ[i, j]

                for (k, l) in get_neighbor_indices(model, (i, j))
                    if (k, l) ∉ cluster
                        neighbor_spin = model.σ[k, l]

                        if rand(model.rng) <
                           1 -
                           exp(2 * model.params.β * dot(spin, r) * dot(neighbor_spin, r))
                            push!(stack, (k, l))
                            push!(cluster, (k, l))
                            flip_spin!(model, (k, l), r)
                        end
                    end
                end
            end

            t > t₀ && (model.observables.avg_cluster_size += length(cluster))
        end

        t > t₀ && update_observables!(model)
    end

    model.observables.avg_cluster_size /= (model.params.n_sweeps * n_clusters)
    compute_observables_statistics!(model)
end

"""
    swendsen_wang!(model::Ising)

Implementation of the Swendsen-Wang multiple cluster algorithm for the Ising model.

Arguments
* `model`: Ising type
"""
function swendsen_wang!(model::Ising)
    P = 1 - exp(-2 * model.params.β)
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)

    for t = 1:model.params.n_sweeps
        clustered = Set{Tuple{Int,Int}}()

        for j = 1:model.params.L, i in model.params.L
            (i, j) ∈ clustered && continue
            old_spin = model.σ[i, j]
            new_spin = rand(model.rng, Int8[-1, 1])
            model.σ[i, j] = new_spin
            push!(clustered, (i, j))

            stack = [(i, j)]
            while length(stack) > 0
                (i, j) = pop!(stack)
                for (k, l) in get_neighbor_indices(model, (i, j))
                    if model.σ[k, l] == old_spin &&
                       (k, l) ∉ clustered &&
                       rand(model.rng) < P
                        model.σ[k, l] = new_spin
                        push!(stack, (k, l))
                        push!(clustered, (k, l))
                    end
                end
            end
        end

        t > t₀ && update_observables!(model)
    end

    compute_observables_statistics!(model)
end

"""
    swendsen_wang!(model::Potts)

"""
function swendsen_wang!(model::Potts)
end

"""
    swendsen_wang!(model::XY)

"""
function swendsen_wang!(model::XY)
end


text:
using IncrementalInference, KernelDensityEstimate

fg = emptyFactorGraph()

N=100

doors = [-100.0;0.0;100.0;300.0]'
pd = kde!(doors,[3.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
doors2 = getPoints(pd);
v1 = addNode!(fg,"x1",doors,N=N)
f1  = addFactor!(fg,[v1],Obsv2( doors2, bws', [1.0]))

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,"x2", tem, N=N)
addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

# monocular sighting would look something like
#addFactor!(fg, Mono, ["x3","l1"], [14.0], [1.0], [1.0])
#addFactor!(fg, Mono, ["x4","l1"], [11.0], [1.0], [1.0])

v3=addNode!(fg,"x3",4.0*randn(1,N)+getVal(v2)+50.0, N=N)
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))
f2 = addFactor!(fg,[v3], Obsv2(doors2, bws', [1.0]))

v4=addNode!(fg,"x4",2.0*randn(1,N)+getVal(v3)+50.0, N=N)
addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

    l1=addNode!(fg, "l1", 0.5*randn(1,N)+getVal(v3)+64.0, N=N)
    addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))
    addFactor!(fg, [v4,l1], Ranged([16.0],[0.5],[1.0]))

v5=addNode!(fg,"x5",2.0*randn(1,N)+getVal(v4)+50.0, N=N)
addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

v6=addNode!(fg,"x6",1.25*randn(1,N)+getVal(v5)+40.0, N=N)
addFactor!(fg,[v5;v6],Odo([40.0]',[1.25]',[1.0]))

v7=addNode!(fg,"x7",2.0*randn(1,N)+getVal(v6) +60.0, N=N)
addFactor!(fg,[v6;v7],Odo([60.0]',[2.0]',[1.0]))

f3 = addFactor!(fg,[v7], Obsv2(doors, bws', [1.0]))

# HMM computed ground truth, extended for 7 poses with landmark
gt = Dict{ASCIIString, Array{Float64,2}}()
gt["x1"]=([0.0;1.97304 ]')' # -0.0342366
gt["x2"]=([50.0; 2.83153 ]')' # 49.8797
gt["x3"]=([100.0; 1.65557 ]')' # 99.8351
gt["x4"]=([150.0; 1.64945 ]')' # 148.637
gt["x5"]=([200.0; 1.77992 ]')' # 198.62
gt["x6"]=([240.0; 2.20466 ]')' # 238.492
gt["x7"]=([300.0; 2.14353 ]')' # 298.467
gt["l1"]=([165.0; 1.17284 ]')' # 164.102

tree = prepBatchTree!(fg);

# list vertices in fg
@show xx,ll = ls(fg)

# do belief propagation inference over tree once
inferOverTreeR!(fg, tree)

inferOverTree!(fg, tree)
println("Inference finished")


text:
mutable struct KSTestResult <: TestResult
    n_obs::Int
    statistic::Float64
    pval::Float64
end

function show(io::IO, kstr::KSTestResult)
    println(io, "\n        One Sample Kolmogorov-Smirnov test")
    println(io, "\n Result:\n", coeftable(kstr))
end

function coeftable(kstr::KSTestResult)
    CoefTable(
        [kstr.n_obs, kstr.statistic, kstr.pval],
        ["N", "D", "P-value"],
        [""],
        3, 2
    )
end

function pkstwo(x::Float64, tol = 1e-06)
    k_max = floor(Int, sqrt(2.0 - log(tol)))
    if x < 1
        z = - (π^2 / 8) / (x^2)
        w = log(x)
        s = 0.0
        k = 1
        while k < k_max
            s += exp(k^2 * z - w)
            k += 2
        end
        return s * sqrt(2 * π)
    else
        z = -2 * x²
        s = -1.0
        old = 0.0
        new = 1.0
        while abs(old - new) > tol
            old = new
            new += 2 * s * exp(z * k^2)
            s *= -1.0
            k += 1
        end
        return new
    end
end

function ccdfNonExact(d::KSDist, x::Float64)
    return 1 - pkstwo(sqrt(d.n) * x)
end

"""
    ks_test(x, distr)

One sample Kolmogorov-Smirnov test

The argument `x` can be a numeric `Vector`.

The argument `distr` can be a `Distributions.UnivariateDistribution`.
"""
function ks_test(
    x::NumVector,
    distr::UnivariateDistribution;
    exact::Union{Bool, Nothing}=nothing
)::KSTestResult
    n = length(x)
    ties = false
    if length(unique(x)) < n
        @warn "ties should not be present for the Kolmogorov-Smirnov test"
        ties = true
    end

    if isnothing(exact)
        exact = (n < 100) & !ties
    end

    x_density_sorted = D.cdf.(distr, sort(x)) - collect(0:(n-1))/n
    statistic = max(
        Base.maximum(x_density_sorted),
        Base.maximum(1/n .- x_density_sorted)
    )

    if exact
        pval = D.ccdf(KSDist(n), statistic)
    else
        pval = ccdfNonExact(KSDist(n), statistic)
    end

    res = KSTestResult(n, statistic, pval)
    return res
end


text:
##########
##########
# Definition of functions used to look at species rewiring in relevant paper
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-04-29
##########
##########

# set up =======================================================================

using LightGraphs
using PyCall

run(`$(PyCall.python) -m pip install python-igraph`)
run(`$(PyCall.python) -m pip install partition_igraph`)

ig = pyimport("igraph");
pyimport("partition_igraph");
# main function ================================================================

"""
    network_metrics(species_list, network)

Get all metrics of interest from the focal network and put them in a vector
"""
function network_metrics(species_list::Array, network::Array)
    # simple metrics

    # initialize vector of metrics for each species
    species_network_metrics = Vector{Any}(undef, length(species_list[:,1]))

    # order
    S = length(species_list[:,1])
    # interactions
    L = sum(network)
    # linkage density
    SL = L/S
    # Connectance
    C = L/(S*(S-1))
    # degree
    all_in_deg = Vector{Any}(undef, length(species_list[:,1]))
    all_out_deg = Vector{Any}(undef, length(species_list[:,1]))
    all_total_deg = Vector{Any}(undef, length(species_list[:,1]))
    for i in 1:length(species_list[:,1])
        in_deg = sum(network[i,:])
        out_deg = sum(network[:,i])
        total_deg = in_deg + out_deg
        species_network_metrics[i] = [Float64(i),
                                      Float64(in_deg),
                                      Float64(out_deg),
                                      Float64(total_deg)]
        all_in_deg[i] = in_deg
        all_out_deg[i] = out_deg
        all_total_deg[i] = total_deg
    end
    unique_total_deg = unique(all_total_deg)
    # degree distributions
    unique_in_deg = unique(all_in_deg)
    unique_out_deg = unique(all_out_deg)
    pk_in = Vector{Any}(undef, length(unique_in_deg))
    pk_out = Vector{Any}(undef, length(unique_out_deg))
    j = 1
    for i in unique_in_deg # assuming each i is each k
        pk_in_temp = (length(all_in_deg[all_in_deg .== i]))/S # P(k) = N(k)/S
        pk_in[j] = [i, pk_in_temp]
        j += 1
    end
    j = 1
    for i in unique_out_deg # assuming each i is each k
        pk_out_temp = (length(all_out_deg[all_out_deg .== i]))/S # P(k) = N(k)/S
        pk_out[j] = [i, pk_out_temp]
        j += 1
    end
    # diameter
    interaction_matrix = convert(Array{Int64}, network)
    graph = DiGraph(interaction_matrix)
    diameter = LightGraphs.diameter(SimpleGraph(graph))
    # local clustering coefficient
    for i in 1:length(species_network_metrics)
        cc_i = local_clustering_coefficient(graph, i)
        append!(species_network_metrics[i], cc_i)
    end
    # local (by degree) clustering coefficient
    cc_k = Vector{Any}(undef, length(unique_total_deg))
    for k in 1:length(unique_total_deg)
        cc_1 = []
        for i in 1:length(species_network_metrics)
            if species_network_metrics[i][4] == unique_total_deg[k] # position 4 is total deg
                append!(cc_1, species_network_metrics[i][5])
            end
        end
        cc_k[k] = [unique_total_deg[k], (sum(cc_1)/S)]
    end
    # global clustering coefficient
    global_cc = global_clustering_coefficient(graph)
    # modularity
    #simple_graph = SimpleGraph(graph)
    #modularity(simple_graph)
    # nestedness

    # web metrics
    web_metrics = [S, L, SL, C, diameter, global_cc]

    return web_metrics, all_in_deg, all_out_deg, all_total_deg

end

## practice with the example:
# https://www.juliabloggers.com/a-tutorial-on-igraph-for-julia/

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
#====================================================================================================
This script generates latex symbols and names from the Julia RPL list and outputs a lua file
Based on https://gitlab.com/ExpandingMan/cmp-latex/-/blob/master/generate_items.jl
====================================================================================================#

using REPL

function add!(s, pair)
    x, y = pair
    if count(elem -> elem[1] == x, s) == 0
        push!(s, string(x) => string(y))
    end
end

function unicode_list()
    s = []
    out = readlines(open(`curl -L http://milde.users.sourceforge.net/LUCR/Math/data/unimathsymbols.txt`))
    out = filter(x -> !startswith(x, "#"), out) # Remove comments
    out = split.(out, '^') # Split on ^
    out = filter( x -> startswith(x[3], '\\') || startswith(x[4], '\\'), out) # Only get some elements
    for elem in out
        if strip(elem[2]) == ""
            continue
        end
        if startswith(elem[3], "\\")
            add!(s, elem[3] => elem[2])
        end
        if startswith(elem[4], "\\")
            add!(s, elem[4] => elem[2])
        end
        m = match(r"= (\\\w+)", elem[end])
        if !isnothing(m)
            add!(s, m.captures[1] => elem[2])
        end
        m = match(r"# (\\\w+)", elem[end])
        if !isnothing(m)
            add!(s, m.captures[1] => elem[2])
        end
    end
    _s = collect(REPL.REPLCompletions.latex_symbols)
    for elem in _s
        add!(s, elem)
    end
    add!(s, "\\mathord" => "⍹")
    return s
end

function luaitem(io::IO, p::Pair)
    k, v = p
    k = replace(k, "\\"=>"\\\\")
    kq = "\"$k\""
    vq = "\"$v\""
    label = '"'*k*" "*v*'"'
    write(io, "{word=$kq, label=$label, insertText=$vq, filterText=$kq}")
end

function luaitems(io::IO, l=unicode_list())
    foreach(l) do p
        print(io, "  ")
        luaitem(io, p)
        print(io, ",\n")
    end
end

function luafile(io::IO, l=unicode_list(), objname="symbols")
    print(io, "local $objname = {\n")
    luaitems(io, l)
    print(io, "}\n")
    print(io, "\n")
    print(io, "return $objname")
end

luafile(fname::AbstractString="./lua/cmp_latex_symbols/items.lua", l=unicode_list()) = open(io -> luafile(io, l), fname, write=true, create=true)

max_stars_repo_name: JC-Doyle/cmp-latex-symbols

max_issues_repo_name: JC-Doyle/cmp-latex-symbols

max_forks_repo_name: JC-Doyle/cmp-latex-symbols


text:
using CompScienceMeshes
using BEAST

Faces = meshsphere(1.0, 0.35)
Edges = skeleton(Faces,1)

faces = barycentric_refinement(Faces)
edges = skeleton(faces,1)
# verts = skeleton(faces,0)

E = 1
Edge = cells(Edges)[1]

port_idx = numvertices(Faces) + E
ptch_idx = Edge[1]

# patch = Mesh(vertices(faces), filter(c -> ptch_idx in c, cells(faces)))
patch_idcs = Int[]
for (i,face) in enumerate(cells(faces))
    if ptch_idx in face
        push!(patch_idcs, i)
    end
end
patch = Mesh(vertices(faces), cells(faces)[patch_idcs])

port = Mesh(vertices(edges), filter(c -> port_idx in c, cells(boundary(patch))))

@show numcells(patch)
@show numcells(port)

# D, C, d, c, d0, d1,
RT_int, RT_prt, x_int, x_prt = BEAST.buildhalfbc2(patch, port, nothing)

BF = BEAST.Shape{Float64}[]
for (m,bf) in enumerate(RT_int.fns)
    for sh in bf
        cellid = patch_idcs[sh.cellid]
        BEAST.add!(BF,cellid, sh.refid, sh.coeff * x_int[m])
    end
end

for (m,bf) in enumerate(RT_prt.fns)
    for sh in bf
        cellid = patch_idcs[sh.cellid]
        BEAST.add!(BF,cellid, sh.refid, sh.coeff * x_prt[m])
    end
end

# RT_prt = raviartthomas(patch, cellpairs(patch, port))
# Lx = lagrangecxd0(patch)
# @assert numfunctions(Lx) == numcells(patch)
#
# Id = BEAST.Identity()
# div_RT_prt = divergence(RT_prt)
# Z, store = BEAST.allocatestorage(Id, Lx, div_RT_prt, Val{:bandedstorage}, BEAST.LongDelays{:ignore})
# BEAST.assemble_local_mixed!(Id, Lx, div_RT_prt, store)

bcs3 = BEAST.buffachristiansen3(Faces)
bcs2 = BEAST.buffachristiansen2(Faces)
bcs1 = BEAST.buffachristiansen(Faces)
rts = BEAST.raviartthomas(Faces)

G1 = assemble(BEAST.NCross(), bcs1, rts)
Q1 = assemble(BEAST.Identity(), divergence(bcs1), divergence(bcs1))

G2 = assemble(BEAST.NCross(), bcs2, rts)
Q2 = assemble(BEAST.Identity(), divergence(bcs2), divergence(bcs2))

G3 = assemble(BEAST.NCross(), bcs3, rts)
Q3 = assemble(BEAST.Identity(), divergence(bcs3), divergence(bcs3))

using LinearAlgebra
@show cond(G3)

using LinearAlgebra
@assert (numcells(Faces) - 1) == (size(Q,1) - rank(Q))
@assert cond(G) < 3.5

function compress!(space)
    T = scalartype(space)
    for (i,fn) in pairs(space.fns)
        shapes = Dict{Tuple{Int,Int},T}()
        for shape in fn
            v = get(shapes, (shape.cellid, shape.refid), zero(T))
            shapes[(shape.cellid, shape.refid)] = v + shape.coeff
            # set!(shapes, (shape.cellid, shape.refid), v + shape.coeff)
        end
        space.fns[i] = [BEAST.Shape(k[1], k[2], v) for (k,v) in shapes]
    end
end

getch(geo,i) = chart(geo, cells(geo)[i])

import PlotlyJS
function showfn(space,i)
    geo = geometry(space)
    T = coordtype(geo)
    X = Dict{Int,T}()
    Y = Dict{Int,T}()
    Z = Dict{Int,T}()
    U = Dict{Int,T}()
    V = Dict{Int,T}()
    W = Dict{Int,T}()
    for sh in space.fns[i]
        chrt = chart(geo, cells(geo)[sh.cellid])
        nbd = center(chrt)
        vals = refspace(space)(nbd)
        x,y,z = cartesian(nbd)
        @show vals[sh.refid].value
        u,v,w = vals[sh.refid].value
        # @show x, y, z
        # @show u, v, w
        X[sh.cellid] = x
        Y[sh.cellid] = y
        Z[sh.cellid] = z
        U[sh.cellid] = get(U,sh.cellid,zero(T)) + sh.coeff * u
        V[sh.cellid] = get(V,sh.cellid,zero(T)) + sh.coeff * v
        W[sh.cellid] = get(W,sh.cellid,zero(T)) + sh.coeff * w
    end
    X = collect(values(X))
    Y = collect(values(Y))
    Z = collect(values(Z))
    U = collect(values(U))
    V = collect(values(V))
    W = collect(values(W))
    PlotlyJS.cone(x=X,y=Y,z=Z,u=U,v=V,w=W)
end