text:
collect(enumerate(VE)) = [
 (v1, [e1, e4, e14]),
 (v2, [e2, e5, e15]),
 (v3, [e1, e2, e7, e17]),
 (v4, [e3, e4, e19]),
 (v5, [e3, e6, e20]),
 (v6, [e5, e6, e8, e21]),
 (v7, [e7, e9, e46, e47]),
 (v8, [e8, e9, e54, e55]),
 (v9, [e10, e12, e14]),
 (v10, [e10, e13, e16]),
 (v11, [e11, e12, e19]),
 (v12, [e11, e13, e20]),
 (v13, [e15, e16, e18, e23]),
 (v14, [e17, e18, e48, e49]),
 (v15, [e21, e22, e53, e54]),
 (v16, [e22, e23, e62, e63]),
 (v17, [e24, e26, e36]),
 (v18, [e24, e27, e37]),
 (v19, [e25, e26, e38]),
 (v20, [e25, e27, e39]),
 (v21, [e28, e30, e36]),
 (v22, [e28, e31, e37]),
 (v23, [e29, e30, e38]),
 (v24, [e29, e33, e39]),
 (v25, [e32, e33, e35, e41]),
 (v26, [e31, e32, e34, e43, e65, e66, e67]),
 (v27, [e34, e35, e60, e61]),
 (v28, [e40, e42, e50, e51, e67]),
 (v29, [e40, e43, e57, e58]),
 (v30, [e41, e42, e44, e45]),
 (v31, [e44, e48, e59]),
 (v32, [e45, e50, e60]),
 (v33, [e46, e49, e62]),
 (v34, [e47, e51, e64]),
 (v35, [e52, e56, e59]),
 (v36, [e52, e57, e61, e65]),
 (v37, [e53, e56, e63]),
 (v38, [e55, e58, e64, e66])]

collect(enumerate(EF)) = [
 (e1, [f1, f4]),
 (e2, [f2, f5]),
 (e3, [f1, f6]),
 (e4, [f1, f7]),
 (e5, [f2, f8]),
 (e6, [f1, f9]),
 (e7, [f1, f2, f18, f20]),
 (e8, [f1, f2, f21, f22]),
 (e9, [f1, f2, f25, f26]),
 (e10, [f3, f4]),
 (e11, [f3, f6]),
 (e12, [f3, f7]),
 (e13, [f3, f9]),
 (e14, [f4, f7]),
 (e15, [f5, f8]),
 (e16, [f4, f9]),
 (e17, [f4, f5, f18, f20]),
 (e18, [f4, f5, f28, f29]),
 (e19, [f6, f7]),
 (e20, [f6, f9]),
 (e21, [f8, f9, f21, f22]),
 (e22, [f8, f9, f25, f27]),
 (e23, [f8, f9, f28, f29]),
 (e24, [f10, f13]),
 (e25, [f10, f14]),
 (e26, [f10, f15]),
 (e27, [f10, f16]),
 (e28, [f11, f13]),
 (e29, [f11, f14]),
 (e30, [f11, f15]),
 (e31, [f11, f16]),
 (e32, [f12, f17]),
 (e33, [f11, f16]),
 (e34, [f11, f12, f30, f32]),
 (e35, [f11, f12, f23, f24]),
 (e36, [f13, f15]),
 (e37, [f13, f16]),
 (e38, [f14, f15]),
 (e39, [f14, f16]),
 (e40, [f16, f17]),
 (e41, [f16, f17, f23, f24]),
 (e42, [f16, f17, f18, f19]),
 (e43, [f16, f17]),
 (e44, [f18, f23]),
 (e45, [f19, f24]),
 (e46, [f20, f25]),
 (e47, [f18, f26]),
 (e48, [f18, f28]),
 (e49, [f20, f29]),
 (e50, [f19, f30]),
 (e51, [f18, f31]),
 (e52, [f21, f23]),
 (e53, [f21, f27]),
 (e54, [f22, f25]),
 (e55, [f21, f26]),
 (e56, [f21, f28]),
 (e57, [f21]),
 (e58, [f21]),
 (e59, [f23, f28]),
 (e60, [f24, f30]),
 (e61, [f23, f32]),
 (e62, [f25, f29]),
 (e63, [f27, f28]),
 (e64, [f26, f31]),
 (e65, [f32]),
 (e66, [f31]),
 (e67, [f30, f31])]

collect(enumerate(EV)) = [
 (e1, [v1, v3]),
 (e2, [v2, v3]),
 (e3, [v4, v5]),
 (e4, [v1, v4]),
 (e5, [v2, v6]),
 (e6, [v5, v6]),
 (e7, [v3, v7]),
 (e8, [v6, v8]),
 (e9, [v7, v8]),
 (e10, [v9, v10]),
 (e11, [v11, v12]),
 (e12, [v9, v11]),
 (e13, [v10, v12]),
 (e14, [v1, v9]),
 (e15, [v2, v13]),
 (e16, [v10, v13]),
 (e17, [v3, v14]),
 (e18, [v13, v14]),
 (e19, [v4, v11]),
 (e20, [v5, v12]),
 (e21, [v6, v15]),
 (e22, [v15, v16]),
 (e23, [v13, v16]),
 (e24, [v17, v18]),
 (e25, [v19, v20]),
 (e26, [v17, v19]),
 (e27, [v18, v20]),
 (e28, [v21, v22]),
 (e29, [v23, v24]),
 (e30, [v21, v23]),
 (e31, [v22, v26]),
 (e32, [v25, v26]),
 (e33, [v24, v25]),
 (e34, [v26, v27]),
 (e35, [v25, v27]),
 (e36, [v17, v21]),
 (e37, [v18, v22]),
 (e38, [v19, v23]),
 (e39, [v20, v24]),
 (e40, [v28, v29]),
 (e41, [v25, v30]),
 (e42, [v28, v30]),
 (e43, [v26, v29]),
 (e44, [v30, v31]),
 (e45, [v30, v32]),
 (e46, [v7, v33]),
 (e47, [v7, v34]),
 (e48, [v14, v31]),
 (e49, [v14, v33]),
 (e50, [v28, v32]),
 (e51, [v28, v34]),
 (e52, [v35, v36]),
 (e53, [v15, v37]),
 (e54, [v8, v15]),
 (e55, [v8, v38]),
 (e56, [v35, v37]),
 (e57, [v29, v36]),
 (e58, [v29, v38]),
 (e59, [v31, v35]),
 (e60, [v27, v32]),
 (e61, [v27, v36]),
 (e62, [v16, v33]),
 (e63, [v16, v37]),
 (e64, [v34, v38]),
 (e65, [v26, v36]),
 (e66, [v26, v38]),
 (e67, [v26, v28])]

 collect(enumerate(FV)) = [
 (f1, [v7, v4, v3, v5, v8, v6, v1]),
 (f2, [v7, v2, v3, v8, v6]),
 (f3, [v9, v10, v11, v12]),
 (f4, [v9, v13, v10, v14, v3, v1]),
 (f5, [v13, v14, v2, v3]),
 (f6, [v4, v11, v5, v12]),
 (f7, [v4, v9, v11, v1]),
 (f8, [v13, v2, v16, v15, v6]),
 (f9, [v13, v10, v16, v5, v15, v6, v12]),
 (f10, [v20, v19, v17, v18]),
 (f11, [v23, v26, v25, v27, v21, v22, v24]),
 (f12, [v26, v25, v27]),
 (f13, [v17, v21, v22, v18]),
 (f14, [v23, v24, v19, v20]),
 (f15, [v23, v19, v17, v21]),
 (f16, [v24, v30, v28, v20, v26, v25, v29, v22, v18]),
 (f17, [v26, v25, v29, v28, v30]),
 (f18, [v7, v14, v31, v34, v3, v28, v30]),
 (f19, [v32, v28, v30]),
 (f20, [v7, v14, v3, v33]),
 (f21, [v37, v35, v29, v36, v8, v15, v6, v38]),
 (f22, [v8, v15, v6]),
 (f23, [v25, v31, v35, v27, v36, v30]),
 (f24, [v25, v32, v27, v30]),
 (f25, [v7, v16, v8, v15, v33]),
 (f26, [v7, v34, v8, v38]),
 (f27, [v37, v16, v15]),
 (f28, [v37, v13, v14, v31, v35, v16]),
 (f29, [v13, v14, v16, v33]),
 (f30, [v26, v32, v27, v28]),
 (f31, [v26, v34, v28, v38]),
 (f32, [v26, v27, v36])]

collect(enumerate(Lar.cop2lar(copFE))) = [
 (f1, [e1, e3, e4, e6, e7, e8, e9]),
 (f2, [e2, e5, e7, e8, e9]),
 (f3, [e10, e11, e12, e13]),
 (f4, [e1, e10, e14, e16, e17, e18]),
 (f5, [e2, e15, e17, e18]),
 (f6, [e3, e11, e19, e20]),
 (f7, [e4, e12, e14, e19]),
 (f8, [e5, e15, e21, e22, e23]),
 (f9, [e6, e13, e16, e20, e21, e22, e23]),
 (f10, [e24, e25, e26, e27]),
 (f11, [e28, e29, e30, e31, e33, e34, e35]),
 (f12, [e32, e34, e35]),
 (f13, [e24, e28, e36, e37]),
 (f14, [e25, e29, e38, e39]),
 (f15, [e26, e30, e36, e38]),
 (f16, [e27, e31, e33, e37, e39, e40, e41, e42, e43]),
 (f17, [e32, e40, e41, e42, e43]),
 (f18, [e7, e17, e42, e44, e47, e48, e51]),
 (f19, [e42, e45, e50]),
 (f20, [e7, e17, e46, e49]),
 (f21, [e8, e21, e52, e53, e55, e56, e57, e58]),
 (f22, [e8, e21, e54]),
 (f23, [e35, e41, e44, e52, e59, e61]),
 (f24, [e35, e41, e45, e60]),
 (f25, [e9, e22, e46, e54, e62]),
 (f26, [e9, e47, e55, e64]),
 (f27, [e22, e53, e63]),
 (f28, [e18, e23, e48, e56, e59, e63]),
 (f29, [e18, e23, e49, e62]),
 (f30, [e34, e50, e60, e67]),
 (f31, [e51, e64, e66, e67]),
 (f32, [e34, e61, e65])]

FF = copFV * copFV'
ff = [(SparseArrays.findnz(FF[k, :]))[1] for k = 1:size(FF, 2)];
collect(enumerate(ff)) = [
 (1, [1, 2, 4, 5, 6, 7, 8, 9, 18, 20, 21, 22, 25, 26]),
 (2, [1, 2, 4, 5, 8, 9, 18, 20, 21, 22, 25, 26]),
 (3, [3, 4, 6, 7, 9]),
 (4, [1, 2, 3, 4, 5, 7, 8, 9, 18, 20, 28, 29]),
 (5, [1, 2, 4, 5, 8, 9, 18, 20, 28, 29]),
 (6, [1, 3, 6, 7, 9]),
 (7, [1, 3, 4, 6, 7]),
 (8, [1, 2, 4, 5, 8, 9, 21, 22, 25, 27, 28, 29]),
 (9, [1, 2, 3, 4, 5, 6, 8, 9, 21, 22, 25, 27, 28, 29]),
 (10, [10, 13, 14, 15, 16]),
 (11, [11, 12, 13, 14, 15, 16, 17, 23, 24, 30, 31, 32]),
 (12, [11, 12, 16, 17, 23, 24, 30, 31, 32]),
 (13, [10, 11, 13, 15, 16]),
 (14, [10, 11, 14, 15, 16]),
 (15, [10, 11, 13, 14, 15]),
 (16, [10, 11, 12, 13, 14, 16, 17, 18, 19, 21, 23, 24, 30, 31, 32]),
 (17, [11, 12, 16, 17, 18, 19, 21, 23, 24, 30, 31, 32]),
 (18, [1, 2, 4, 5, 16, 17, 18, 19, 20, 23, 24, 25, 26, 28, 29, 30, 31]),
 (19, [16, 17, 18, 19, 23, 24, 30, 31]),
 (20, [1, 2, 4, 5, 18, 20, 25, 26, 28, 29]),
 (21, [1, 2, 8, 9, 16, 17, 21, 22, 23, 25, 26, 27, 28, 31, 32]),
 (22, [1, 2, 8, 9, 21, 22, 25, 26, 27]),
 (23, [11, 12, 16, 17, 18, 19, 21, 23, 24, 28, 30, 32]),
 (24, [11, 12, 16, 17, 18, 19, 23, 24, 30, 32]),
 (25, [1, 2, 8, 9, 18, 20, 21, 22, 25, 26, 27, 28, 29]),
 (26, [1, 2, 18, 20, 21, 22, 25, 26, 31]),
 (27, [8, 9, 21, 22, 25, 27, 28, 29]),
 (28, [4, 5, 8, 9, 18, 20, 21, 23, 25, 27, 28, 29]),
 (29, [4, 5, 8, 9, 18, 20, 25, 27, 28, 29]),
 (30, [11, 12, 16, 17, 18, 19, 23, 24, 30, 31, 32]),
 (31, [11, 12, 16, 17, 18, 19, 21, 26, 30, 31, 32]),
 (32, [11, 12, 16, 17, 21, 23, 24, 30, 31, 32])]

for (V,EV) in faces2D
	VV = [[k] for k=1:size(V,2)]
	model = (convert(Lar.Points,V), [VV,EV,FV])::Lar.LARmodel
	println(EV)
	GL.VIEW(GL.numbering(.1)(model, GL.COLORS[1], 0.1));
end

using ViewerGL, LinearAlgebra, SparseArrays, LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation; GL = ViewerGL

V = [0.6540617999999998 0.2054992 0.2972308000000001; 0.7142364999999998 0.1455625 0.969203; 0.6872816200343089 0.1724107693563844 0.6681972618528844; 0.5941250999999997 0.8769965 0.36249240000000005; 0.6542998079614951 0.8170598007732168 1.0344646993560203; 0.6929229746417847 0.3843473621810608 0.9924101013070119; 0.6660927570227951 0.3906268097799347 0.66824889023436; 0.6924256885846046 0.3894687221287142 0.9924113129852241; 1.3260341 0.2707609000000001 0.2428771; 1.3862088 0.21082410000000004 0.9148494000000003; 1.2660973 0.9422582000000003 0.3081388; 1.3262719920385087 0.8823213992267827 0.9801110006439793; 0.7946185131077457 0.1533691699791747 0.9627011677876783; 0.795311803386693 0.18207414025178234 0.6679835359595829; 0.6978893364451173 0.3892254686183435 0.9924004714668809; 0.8049461769791598 0.3844586322271989 0.9823884300888783; -0.38740630000000004 0.49022260000000006 0.45363390000000015; 0.3249123 0.707347 0.5231232000000001; -0.17028190000000007 -0.08642419999999995 0.029717700000000125; 0.542036694733776 0.13070013212634024 0.09920695360173326; -0.317917 0.06630639999999999 1.0658723; 0.3944015999999999 0.28343080000000004 1.1353616; -0.10079260000000007 -0.5103404 0.6419561000000001; 0.6115259947337761 -0.2932160678736597 0.7114453536017331; 0.4729938678063347 0.07470259548304928 0.9819171231357612; 0.4673971777131105 0.08956648867512645 0.992844167328483; 0.4666328572520818 0.07498584149882567 0.9828352213039846; 0.4707836778272794 0.18222015069807973 0.6686268661789595; 0.4673140236772549 0.08983923498806379 0.9928443441979011; 0.5039882186318663 0.07330882745780876 0.6685562398308738; 0.7899025999999999 0.060579300000000016 0.6679888999999999; 0.46601 0.0749997 0.6686315999999999; 0.8043229999999999 0.38447249999999994 0.6679745999999999; 0.4804303999999999 0.3988929 0.6686172999999999; 0.7905452000000002 0.060564999999999994 0.9922023000000001; 0.46665270000000003 0.0749854 0.9928450000000001; 0.8049656 0.3844582 0.9921880000000001; 0.4810731000000001 0.3988786 0.9928307000000001];
copEV = SparseArrays.sparse([1, 4, 14, 2, 5, 15, 1, 2, 7, 17, 3, 4, 19, 3, 6, 20, 5, 6, 8, 21, 7, 9, 46, 47, 8, 9, 54, 55, 10, 12, 14, 10, 13, 16, 11, 12, 19, 11, 13, 20, 15, 16, 18, 23, 17, 18, 48, 49, 21, 22, 53, 54, 22, 23, 62, 63, 24, 26, 36, 24, 27, 37, 25, 26, 38, 25, 27, 39, 28, 30, 36, 28, 31, 37, 29, 30, 38, 29, 33, 39, 32, 33, 35, 41, 31, 32, 34, 43, 65, 66, 67, 34, 35, 60, 61, 40, 42, 50, 51, 67, 40, 43, 57, 58, 41, 42, 44, 45, 44, 48, 59, 45, 50, 60, 46, 49, 62, 47, 51, 64, 52, 56, 59, 52, 57, 61, 65, 53, 56, 63, 55, 58, 64, 66], [1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 38, 38, 38, 38], Int8[-1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]);
copFE = SparseArrays.sparse([1, 4, 2, 5, 1, 6, 1, 7, 2, 8, 1, 9, 1, 2, 18, 20, 1, 2, 21, 22, 1, 2, 25, 26, 3, 4, 3, 6, 3, 7, 3, 9, 4, 7, 5, 8, 4, 9, 4, 5, 18, 20, 4, 5, 28, 29, 6, 7, 6, 9, 8, 9, 21, 22, 8, 9, 25, 27, 8, 9, 28, 29, 10, 13, 10, 14, 10, 15, 10, 16, 11, 13, 11, 14, 11, 15, 11, 16, 12, 17, 11, 16, 11, 12, 30, 32, 11, 12, 23, 24, 13, 15, 13, 16, 14, 15, 14, 16, 16, 17, 16, 17, 23, 24, 16, 17, 18, 19, 16, 17, 18, 23, 19, 24, 20, 25, 18, 26, 18, 28, 20, 29, 19, 30, 18, 31, 21, 23, 21, 27, 22, 25, 21, 26, 21, 28, 21, 21, 23, 28, 24, 30, 23, 32, 25, 29, 27, 28, 26, 31, 32, 31, 30, 31], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22, 23, 23, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 37, 37, 38, 38, 39, 39, 40, 40, 41, 41, 41, 41, 42, 42, 42, 42, 43, 43, 44, 44, 45, 45, 46, 46, 47, 47, 48, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54, 54, 55, 55, 56, 56, 57, 58, 59, 59, 60, 60, 61, 61, 62, 62, 63, 63, 64, 64, 65, 66, 67, 67], Int8[1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1]);
size(V), size(copEV), size(copFE)

using SparseArrays
V = [0.6540617999999998 0.7142364999999998 0.6872816200343089 0.5941250999999997 0.6542998079614951 0.6929229746417847 0.6660927570227951 0.6924256885846046 1.3260341 1.3862088 1.2660973 1.3262719920385087 0.7946185131077457 0.795311803386693 0.8049461769791598 0.6978893364451173 -0.38740630000000004 0.3249123 -0.17028190000000007 0.542036694733776 -0.317917 0.3944015999999999 -0.10079260000000007 0.6115259947337761 0.4729938678063347 0.4673971777131105 0.4666328572520818 0.5039882186318663 0.4707836778272794 0.4673140236772549 0.7899025999999999 0.46601 0.8043229999999999 0.4804303999999999 0.7905452000000002 0.46665270000000003 0.8049656 0.4810731000000001; 0.2054992 0.1455625 0.1724107693563844 0.8769965 0.8170598007732168 0.3843473621810608 0.3906268097799347 0.3894687221287142 0.2707609000000001 0.21082410000000004 0.9422582000000003 0.8823213992267827 0.1533691699791747 0.18207414025178234 0.3844586322271989 0.3892254686183435 0.49022260000000006 0.707347 -0.08642419999999995 0.13070013212634024 0.06630639999999999 0.28343080000000004 -0.5103404 -0.2932160678736597 0.07470259548304928 0.08956648867512645 0.07498584149882567 0.07330882745780876 0.18222015069807973 0.08983923498806379 0.060579300000000016 0.0749997 0.38447249999999994 0.3988929 0.060564999999999994 0.0749854 0.3844582 0.3988786; 0.2972308000000001 0.969203 0.6681972618528844 0.36249240000000005 1.0344646993560203 0.9924101013070119 0.66824889023436 0.9924113129852241 0.2428771 0.9148494000000003 0.3081388 0.9801110006439793 0.9627011677876783 0.6679835359595829 0.9823884300888783 0.9924004714668809 0.45363390000000015 0.5231232000000001 0.029717700000000125 0.09920695360173326 1.0658723 1.1353616 0.6419561000000001 0.7114453536017331 0.9819171231357612 0.992844167328483 0.9828352213039846 0.6685562398308738 0.6686268661789595 0.9928443441979011 0.6679888999999999 0.6686315999999999 0.6679745999999999 0.6686172999999999 0.9922023000000001 0.9928450000000001 0.9921880000000001 0.9928307000000001]
VV = [[1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38]]
EV = [[1, 3], [2, 3], [4, 5], [1, 4], [2, 6], [5, 6], [3, 7], [7, 8], [6, 8], [9, 10], [11, 12], [9, 11], [10, 12], [1, 9], [2, 13], [10, 13], [3, 14], [13, 14], [4, 11], [5, 12], [13, 15], [15, 16], [6, 16], [17, 18], [19, 20], [17, 19], [18, 20], [21, 22], [23, 24], [21, 23], [22, 26], [25, 26], [24, 25], [25, 27], [26, 27], [17, 21], [18, 22], [19, 23], [20, 24], [28, 29], [25, 28], [29, 30], [26, 30], [28, 31], [28, 32], [7, 33], [7, 34], [14, 31], [14, 33], [29, 32], [29, 34], [35, 36], [16, 37], [8, 16], [8, 38], [35, 37], [30, 36], [30, 38], [31, 35], [27, 32], [27, 36], [15, 33], [15, 37], [34, 38], [26, 36], [26, 38], [26, 29]]
FV = [[7, 4, 3, 5, 8, 6, 1], [7, 2, 3, 8, 6], [9, 10, 11, 12], [9, 13, 10, 14, 3, 1], [13, 14, 2, 3], [4, 11, 5, 12], [4, 9, 11, 1], [13, 2, 16, 15, 6], [13, 10, 16, 5, 15, 6, 12], [20, 19, 17, 18], [23, 26, 25, 27, 21, 22, 24], [26, 25, 27], [17, 21, 22, 18], [23, 24, 19, 20], [23, 19, 17, 21], [24, 30, 28, 20, 26, 25, 29, 22, 18], [26, 25, 29, 28, 30], [7, 14, 31, 34, 29, 3, 28], [32, 29, 28], [7, 14, 3, 33], [37, 35, 16, 36, 8, 38, 6, 30], [16, 8, 6], [25, 31, 35, 27, 28, 36], [25, 32, 27, 28], [7, 16, 8, 15, 33], [7, 34, 8, 38], [37, 16, 15], [37, 13, 14, 31, 35, 15], [13, 14, 15, 33], [26, 32, 29, 27], [26, 34, 29, 38], [26, 27, 36]]
copFV = Lar.lar2cop(FV)
FF = copFV * copFV'
ff = [(SparseArrays.findnz(FF[k, :]))[1] for k = 1:size(FF, 2)];
ffV = [FV[f] for f in ff]
copFE = SparseArrays.sparse([1, 4, 2, 5, 1, 6, 1, 7, 2, 8, 1, 9, 1, 2, 18, 20, 1, 2, 25, 26, 1, 2, 21, 22, 3, 4, 3, 6, 3, 7, 3, 9, 4, 7, 5, 8, 4, 9, 4, 5, 18, 20, 4, 5, 28, 29, 6, 7, 6, 9, 8, 9, 28, 29, 8, 9, 25, 27, 8, 9, 21, 22, 10, 13, 10, 14, 10, 15, 10, 16, 11, 13, 11, 14, 11, 15, 11, 16, 12, 17, 11, 16, 11, 12, 23, 24, 11, 12, 30, 32, 13, 15, 13, 16, 14, 15, 14, 16, 16, 17, 18, 19, 16, 17, 23, 24, 16, 17, 16, 17, 18, 23, 19, 24, 20, 25, 18, 26, 18, 28, 20, 29, 19, 30, 18, 31, 21, 23, 21, 27, 22, 25, 21, 26, 21, 28, 21, 21, 23, 28, 24, 30, 23, 32, 25, 29, 27, 28, 26, 31, 32, 31, 30, 31], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22, 23, 23, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 37, 37, 38, 38, 39, 39, 40, 40, 40, 40, 41, 41, 41, 41, 42, 42, 43, 43, 44, 44, 45, 45, 46, 46, 47, 47, 48, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54, 54, 55, 55, 56, 56, 57, 58, 59, 59, 60, 60, 61, 61, 62, 62, 63, 63, 64, 64, 65, 66, 67, 67], Int8[1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1]);
FE = Lar.cop2lar(copFE)
# [[1, 3, 4, 6, 7, 8, 9], [2, 5, 7, 8, 9], [10, 11, 12, 13], [1, 10, 14, 16, 17, 18], [2, 15, 17, 18], [3, 11, 19, 20], [4, 12, 14, 19], [5, 15, 21, 22, 23], [6, 13, 16, 20, 21, 22, 23], [24, 25, 26, 27], [28, 29, 30, 31, 33, 34, 35], [32, 34, 35], [24, 28, 36, 37], [25, 29, 38, 39], [26, 30, 36, 38], [27, 31, 33, 37, 39, 40, 41, 42, 43], [32, 40, 41, 42, 43], [7, 17, 40, 44, 47, 48, 51], [40, 45, 50], [7, 17, 46, 49], [9, 23, 52, 53, 55, 56, 57, 58], [9, 23, 54], [34, 41, 44, 52, 59, 61], [34, 41, 45, 60], [8, 22, 46, 54, 62], [8, 47, 55, 64], [22, 53, 63], [18, 21, 48, 56, 59, 63], [18, 21, 49, 62], [35, 50, 60, 67], [51, 64, 66, 67], [35, 61, 65]]
ffE = [sort!(collect(Set(cat(FE[f]...,dims=1)))) for f in ff]
ffE = [FE[f] for f in ff]

for k=1:length(ff)
	fVV = [[v] for v in sort!(collect(Set(cat([FV[f] for f in ff[k]]))))]
	fEV = cat([EV[e] for e in ffE[k]])
	fEV = sort!(Lar.removeDups(Lar.Cells(fEV)))
	fFV = vcat([FV[f] for f in ff[k]])
	#model = ( ([1 0 0.3; 0 1 0.2; 0 0 1] * V)[1:2,:], Lar.Cells[fVV,fEV,fFV]);
	model = (Lar.Points(V'), Lar.Cells[fVV,fEV,fFV]);
	GL.VIEW(GL.numbering(.1)(model, GL.COLORS[1], 0.1));
end


text:
"""
    ℓ = compute_transition_costs(g, t::AbstractArray)
Computes the transition costs `ℓ` given a
function `g` and a time sequence `t`
"""
function compute_transition_costs(g, t::AbstractVector, tol=1e-3)
    T = Float64
    ℓ_lazy = TransitionCostContinuous{T}(g, t)
    N = length(t)
    ℓ = Array{QuadraticForm{T}}(undef, N-1,N)
    for i = 1:N-1
        for ip = i+1:N
            ℓ[i,ip] = ℓ_lazy[i,ip]
        end
    end
    # TODO move/add tests to Lazy, maybe only i to i+1
    for i = 1:N-1, ip = i+1:N
        P, q, r = ℓ[i,ip].P, ℓ[i,ip].q, ℓ[i,ip].r
        minval = -q⋅(P\q)/4+r
        #TODO add epsilon on RHS?
        minval < 0 && error("Transition cost is negative for some values in `compute_transition_costs`, try decreasing rtol.")
    end

    return ℓ
end

"""
    ℓ = compute_transition_costs(g::AbstractArray; t=1:length(g))
Computes the transition costs `ℓ` given a discrete
function `g`.
"""
function compute_discrete_transition_costs(g::AbstractArray, t=1:length(g))
    T = promote_type(eltype(g), Float64)
    ℓ_lazy = TransitionCostDiscrete{T}(g, t)
    N = length(t)
    ℓ = Array{QuadraticForm{T}}(undef, N-1,N)
    for i = 1:N-1
        for ip = i+1:N
            ℓ[i,ip] = ℓ_lazy[i,ip]
        end
    end
    return ℓ
end

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Set tests
isdefined(Main, :OffsetArrays) || @eval Main include("testhelpers/OffsetArrays.jl")
using .Main.OffsetArrays

using Dates

@testset "Construction, collect" begin
    @test Set([1,2,3]) isa Set{Int}
    @test Set{Int}([3]) isa Set{Int}
    data_in = (1,"banana", ())
    s = Set(data_in)
    data_out = collect(s)
    @test s isa Set{Any}
    @test data_out isa Array{Any,1}
    @test all(map(in(data_out), data_in))
    @test length(data_out) == length(data_in)
    let f17741 = x -> x < 0 ? false : 1
        @test isa(Set(x for x = 1:3), Set{Int})
        @test isa(Set(x for x = 1:3 for j = 1:1), Set{Int})
        @test isa(Set(sin(x) for x = 1:3), Set{Float64})
        @test isa(Set(f17741(x) for x = 1:3), Set{Int})
        @test isa(Set(f17741(x) for x = -1:1), Set{Integer})
        @test isa(Set(f17741(x) for x = 1:0), Set{Integer})
    end
    let s1 = Set(["foo", "bar"]), s2 = Set(s1)
        @test s1 == s2
        x = pop!(s1)
        @test s1 != s2
        @test !(x in s1)
        @test x in s2
        push!(s1, "baz")
        push!(s2, "baz2")
        @test "baz" in s1
        @test !("baz" in s2)
        @test !("baz2" in s1)
        @test "baz2" in s2
    end
end

@testset "hash" begin
    s1 = Set(["bar", "foo"])
    s2 = Set(["foo", "bar"])
    s3 = Set(["baz"])
    @test hash(s1) == hash(s2)
    @test hash(s1) != hash(s3)
    d1 = Dict(Set([3]) => 33, Set([2]) => 22)
    d2 = Dict(Set([2]) => 33, Set([3]) => 22)
    @test hash(d1) != hash(d2)
end

@testset "equality" for eq in (isequal, ==)
    @test  eq(Set(), Set())
    @test !eq(Set(), Set([1]))
    @test  eq(Set{Any}(Any[1,2]), Set{Int}([1,2]))
    @test !eq(Set{Any}(Any[1,2]), Set{Int}([1,2,3]))

    # Comparison of unrelated types
    @test  eq(Set{Int}(), Set{AbstractString}())
    @test !eq(Set{Int}(), Set{AbstractString}([""]))
    @test !eq(Set{AbstractString}(), Set{Int}([0]))
    @test !eq(Set{Int}([1]), Set{AbstractString}())
    @test  eq(Set{Any}([1,2,3]), Set{Int}([1,2,3]))
    @test  eq(Set{Int}([1,2,3]), Set{Any}([1,2,3]))
    @test !eq(Set{Any}([1,2,3]), Set{Int}([1,2,3,4]))
    @test !eq(Set{Int}([1,2,3]), Set{Any}([1,2,3,4]))
    @test !eq(Set{Any}([1,2,3,4]), Set{Int}([1,2,3]))
    @test !eq(Set{Int}([1,2,3,4]), Set{Any}([1,2,3]))

    # Special cases
    @test  eq(Set([-0.0]), Set([-0.0]))
    @test !eq(Set([0.0]), Set([-0.0]))
    @test  eq(Set([NaN]), Set([NaN]))
    @test !eq(Set([NaN]), Set([1.0]))
    @test  eq(Set([missing]), Set([missing]))
    @test !eq(Set([missing]), Set([1]))
end

@testset "hash and == for Set/BitSet" begin
    for s = (Set([1]), Set(1:10), Set(-100:7:100))
        b = BitSet(s)
        @test hash(s) == hash(b)
        @test s == b
    end
end

@testset "eltype, empty" begin
    s1 = empty(Set([1,"hello"]))
    @test isequal(s1, Set())
    @test ===(eltype(s1), Any)
    s2 = empty(Set{Float32}([2.0f0,3.0f0,4.0f0]))
    @test isequal(s2, Set())
    @test ===(eltype(s2), Float32)
    s3 = empty(Set([1,"hello"]),Float32)
    @test isequal(s3, Set())
    @test ===(eltype(s3), Float32)
end
@testset "show" begin
    @test sprint(show, Set()) == "Set{Any}()"
    @test repr([Set(),Set()]) == "Set{Any}[Set(), Set()]"
    @test sprint(show, Set(['a'])) == "Set(['a'])"
end
@testset "isempty, length, in, push, pop, delete" begin
    # also test for no duplicates
    s = Set(); push!(s,1); push!(s,2); push!(s,3)
    @test !isempty(s)
    @test in(1,s)
    @test in(2,s)
    @test length(s) == 3
    push!(s,1); push!(s,2); push!(s,3)
    @test length(s) == 3
    @test pop!(s,1) == 1
    @test !in(1,s)
    @test in(2,s)
    @test length(s) == 2
    @test_throws KeyError pop!(s,1)
    @test pop!(s,1,:foo) == :foo
    @test length(delete!(s,2)) == 1
    @test !in(1,s)
    @test !in(2,s)
    @test pop!(s) == 3
    @test length(s) == 0
    @test isempty(s)
    @test_throws ArgumentError pop!(s)
    @test length(Set(['x',120])) == 2
end
@testset "copy" begin
    data_in = (1,2,9,8,4)
    s = Set(data_in)
    c = copy(s)
    @test isequal(s,c)
    v = pop!(s)
    @test !in(v,s)
    @test  in(v,c)
    push!(s,100)
    push!(c,200)
    @test !in(100,c)
    @test !in(200,s)
end

@testset "copy(::KeySet) (issue #41537)" begin
    @test union(keys(Dict(1=>2, 3=>4))) == copy(keys(Dict(1=>2, 3=>4))) == Set([1,3])
end

@testset "copy!" begin
    for S = (Set, BitSet)
        s = S([1, 2])
        for a = ([1], UInt[1], [3, 4, 5], UInt[3, 4, 5])
            @test s === copy!(s, Set(a)) == S(a)
            @test s === copy!(s, BitSet(a)) == S(a)
        end
    end
    s = Set([1, 2])
    s2 = copy(s)
    @test copy!(s, s) == s2
end

@testset "sizehint, empty" begin
    s = Set([1])
    @test isequal(sizehint!(s, 10), Set([1]))
    @test isequal(empty!(s), Set())
    s2 = GenericSet(s)
    sizehint!(s2, 10)
    @test s2 == GenericSet(s)
end
@testset "rehash!" begin
    # Use a pointer type to have defined behavior for uninitialized
    # array element
    s = Set(["a", "b", "c"])
    Base.rehash!(s)
    k = s.dict.keys
    Base.rehash!(s)
    @test length(k) == length(s.dict.keys)
    for i in 1:length(k)
        if isassigned(k, i)
            @test k[i] == s.dict.keys[i]
        else
            @test !isassigned(s.dict.keys, i)
        end
    end
    @test s == Set(["a", "b", "c"])
end

@testset "start, done, next" begin
    for data_in in ((7, 8, 4, 5),
                  ("hello", 23, 2.7, (), [], (1, 8)))
        local data_in, s, t
        s = Set(data_in)

        s_new = Set()
        for el in s
            push!(s_new, el)
        end
        @test isequal(s, s_new)

        t = tuple(s...)
        @test length(t) == length(s)
        for e in t
            @test in(e,s)
        end
    end
end

@testset "union" begin
    for S in (Set, BitSet, Vector)
        s = ∪(S([1,2]), S([3,4]))
        @test s == S([1,2,3,4])
        s = union(S([5,6,7,8]), S([7,8,9]))
        @test s == S([5,6,7,8,9])
        s = S([1,3,5,7])
        union!(s, (2,3,4,5))
        @test s == S([1,3,5,7,2,4]) # order matters for Vector
        let s1 = S([1, 2, 3])
            @test s1 !== union(s1) == s1
            @test s1 !== union(s1, 2:4) == S([1,2,3,4])
            @test s1 !== union(s1, [2,3,4]) == S([1,2,3,4])
            @test s1 !== union(s1, [2,3,4], S([5])) == S([1,2,3,4,5])
            @test s1 === union!(s1, [2,3,4], S([5])) == S([1,2,3,4,5])
        end
    end
    @test union(Set([1]), BitSet()) isa Set{Int}
    @test union(BitSet([1]), Set()) isa BitSet
    @test union([1], BitSet()) isa Vector{Int}
    # union must uniquify
    @test union([1, 2, 1]) == union!([1, 2, 1]) == [1, 2]
    @test union([1, 2, 1], [2, 2]) == union!([1, 2, 1], [2, 2]) == [1, 2]
    s2 = Set([nothing])
    union!(s2, [nothing])
    @test s2 == Set([nothing])

    @testset "promotion" begin
        ints = [1:5, [1, 2], Set([1, 2])]
        floats = [2:0.1:3, [2.0, 3.5], Set([2.0, 3.5])]

        for a in ints, b in floats
            @test eltype(union(a, b)) == Float64
            @test eltype(union(b, a)) == Float64
        end
    end
end

@testset "intersect" begin
    for S in (Set, BitSet, Vector)
        s = S([1,2]) ∩ S([3,4])
        @test s == S()
        s = intersect(S([5,6,7,8]), S([7,8,9]))
        slong = S(collect(3:63))
        # test #36339 length/order short-cut
        @test intersect(S([5,6,7,8]), slong) == intersect(slong, S([5,6,7,8]))
        @test s == S([7,8])
        @test intersect(S([2,3,1]), S([4,2,3]), S([5,4,3,2])) == S([2,3])
        let s1 = S([1,2,3])
            @test s1 !== intersect(s1) == s1
            @test s1 !== intersect(s1, 2:10) == S([2,3])
            @test s1 !== intersect(s1, [2,3,4]) == S([2,3])
            @test s1 !== intersect(s1, [2,3,4], 3:4) == S([3])
            @test s1 === intersect!(s1, [2,3,4], 3:4) == S([3])
        end
    end
    @test intersect(Set([1]), BitSet()) isa Set{Int}
    @test intersect(BitSet([1]), Set()) isa Set{Any}
    @test intersect(BitSet([1]), Set([1])) isa BitSet
    @test intersect(BitSet([1]), Set([1]), Set([1])) isa BitSet
    @test intersect([1], BitSet()) isa Vector{Int}
    # intersect must uniquify
    @test intersect([1, 2, 1]) == intersect!([1, 2, 1]) == [1, 2]
    @test intersect([1, 2, 1], [2, 2]) == intersect!([1, 2, 1], [2, 2]) == [2]

    # issue #25801
    x = () ∩ (:something,)
    y = () ∩ (42,)
    @test isempty(x)
    @test isempty(y)

    # Discussed in PR#41769
    @testset "promotion" begin
        ints = [1:5, [1, 2], Set([1, 2])]
        floats = [2:0.1:3, [2.0, 3.5], Set([2.0, 3.5])]

        for a in ints, b in floats
            @test eltype(intersect(a, b)) == Float64
            @test eltype(intersect(b, a)) == Float64
            @test eltype(intersect(a, a, b)) == Float64
        end
    end

    # 3-argument version is correctly covered
    @test intersect(Set([1,2]), Set([2]), Set([1,2,3])) == Set([2])
    @test intersect(Set([1,2]), Set([2]), Set([1.,2,3])) == Set([2.])
end

@testset "setdiff" begin
    for S in (Set, BitSet, Vector)
        @test setdiff(S([1,2,3]), S())        == S([1,2,3])
        @test setdiff(S([1,2,3]), S([1]))     == S([2,3])
        @test setdiff(S([1,2,3]), S([1,2]))   == S([3])
        @test setdiff(S([1,2,3]), S([1,2,3])) == S()
        @test setdiff(S([1,2,3]), S([4]))     == S([1,2,3])
        @test setdiff(S([1,2,3]), S([4,1]))   == S([2,3])
        let s1 = S([1, 2, 3])
            @test s1 !== setdiff(s1) == s1
            @test s1 !== setdiff(s1, 2:10) == S([1])
            @test s1 !== setdiff(s1, [2,3,4]) == S([1])
            @test s1 !== setdiff(s1, S([2,3,4]), S([1])) == S()
            @test s1 === setdiff!(s1, S([2,3,4]), S([1])) == S()
        end
    end

    @test setdiff(Set([1]), BitSet()) isa Set{Int}
    @test setdiff(BitSet([1]), Set()) isa BitSet
    @test setdiff([1], BitSet()) isa Vector{Int}
    # setdiff must uniquify
    @test setdiff([1, 2, 1]) == setdiff!([1, 2, 1]) == [1, 2]
    @test setdiff([1, 2, 1], [2, 2]) == setdiff!([1, 2, 1], [2, 2]) == [1]

    s = Set([1,3,5,7])
    setdiff!(s,(3,5))
    @test isequal(s,Set([1,7]))
    s = Set([1,2,3,4])
    setdiff!(s, Set([2,4,5,6]))
    @test isequal(s,Set([1,3]))

    # setdiff iterates the shorter set - make sure this algorithm works
    sa, sb = Set([1,2,3,4,5,6,7]), Set([2,3,9])
    @test Set([1,4,5,6,7]) == setdiff(sa, sb) !== sa
    @test Set([1,4,5,6,7]) == setdiff!(sa, sb) === sa
    sa, sb = Set([1,2,3,4,5,6,7]), Set([2,3,9])
    @test Set([9]) == setdiff(sb, sa) !== sb
    @test Set([9]) == setdiff!(sb, sa) === sb
end

@testset "ordering" begin
    @test Set() < Set([1])
    @test Set([1]) < Set([1,2])
    @test !(Set([3]) < Set([1,2]))
    @test !(Set([3]) > Set([1,2]))
    @test Set([1,2,3]) > Set([1,2])
    @test !(Set([3]) <= Set([1,2]))
    @test !(Set([3]) >= Set([1,2]))
    @test Set([1]) <= Set([1,2])
    @test Set([1,2]) <= Set([1,2])
    @test Set([1,2]) >= Set([1,2])
    @test Set([1,2,3]) >= Set([1,2])
    @test !(Set([1,2,3]) >= Set([1,2,4]))
    @test !(Set([1,2,3]) <= Set([1,2,4]))
end

@testset "issubset, symdiff, isdisjoint" begin
    for S in (Set, BitSet, Vector)
        for (l,r) in ((S([1,2]),     S([3,4])),
                      (S([5,6,7,8]), S([7,8,9])),
                      (S([1,2]),     S([3,4])),
                      (S([5,6,7,8]), S([7,8,9])),
                      (S([1,2,3]),   S()),
                      (S(),          S()),
                      (S(),          S([1,2,3])),
                      (S([1,2,3]),   S([1])),
                      (S([1,2,3]),   S([1,2])),
                      (S([1,2,3]),   S([1,2,3])),
                      (S([1,2,3]),   S([4])),
                      (S([1,2,3]),   S([4,1])))
            @test issubset(intersect(l,r), l)
            @test issubset(intersect(l,r), r)
            @test issubset(l, union(l,r))
            @test issubset(r, union(l,r))
            @test isdisjoint(l,l) == isempty(l)
            @test isdisjoint(l,r) == isempty(intersect(l,r))
            if S === Vector
                @test sort(union(intersect(l,r),symdiff(l,r))) == sort(union(l,r))
            else
                @test union(intersect(l,r),symdiff(l,r)) == union(l,r)
            end
        end
        if S !== Vector
            @test ⊆(S([1]), S([1,2]))
            @test ⊊(S([1]), S([1,2]))
            @test !⊊(S([1]), S([1]))
            @test ⊈(S([1]), S([2]))
            @test ⊇(S([1,2]), S([1]))
            @test ⊋(S([1,2]), S([1]))
            @test !⊋(S([1]), S([1]))
            @test ⊉(S([1]), S([2]))
        end
        let s1 = S([1,2,3,4])
            @test s1 !== symdiff(s1) == s1
            @test s1 !== symdiff(s1, S([2,4,5,6])) == S([1,3,5,6])
            @test s1 !== symdiff(s1, S([2,4,5,6]), [1,6,7]) == S([3,5,7])
            @test s1 === symdiff!(s1, S([2,4,5,6]), [1,6,7]) == S([3,5,7])
        end
    end
    @test symdiff(Set([1,2,3,4]), Set([2,4,5,6])) == Set([1,3,5,6])
    @test symdiff(Set([1]), BitSet()) isa Set{Int}
    @test symdiff(BitSet([1]), Set{Int}()) isa BitSet
    @test symdiff([1], BitSet()) isa Vector{Int}
    # symdiff must NOT uniquify
    @test symdiff([1, 2, 1]) == symdiff!([1, 2, 1]) == [2]
    @test symdiff([1, 2, 1], [2, 2]) == symdiff!([1, 2, 1], [2, 2]) == [2]

    # Base.hasfastin
    @test all(Base.hasfastin, Any[Dict(1=>2), Set(1), BitSet(1), 1:9, 1:2:9,
                                  Dict, Set, BitSet, UnitRange, StepRange])
    @test !any(Base.hasfastin, Any[[1, 2, 3], "123",
                                   Array, String])

    # tests for Dict
    d1 = Dict(1=>nothing, 2=>nothing)
    d2 = Dict(1=>nothing, 3=>nothing)
    d3 = Dict(1=>nothing, 2=>nothing, 3=>nothing)
    @test d3 == merge(d1, d2)
    @test !issubset(d1, d2)
    @test !issubset(d2, d1)
    @test !issubset(d3, d1)
    @test !issubset(d3, d2)
    @test issubset(d1, d3)
    @test issubset(d2, d3)

    # no fast in, long enough container
    @test issubset(Set(Bool[]), rand(Bool, 100)) == true
    # neither has a fast in, right doesn't have a length
    @test isdisjoint([1, 3, 5, 7, 9], Iterators.filter(iseven, 1:10))
end

@testset "unique" begin
    u = @inferred(unique([1, 1, 2]))
    @test in(1, u)
    @test in(2, u)
    @test length(u) == 2
    @test unique(iseven, []) == []
    # type promotion
    @test unique(x -> x^2, [1, 3.]) == [1, 3.]
    @test @inferred(unique(iseven, [5, 1, 8, 9, 3, 4, 10, 7, 2, 6])) == [5, 8]
    @test @inferred(unique(x->x^2, Integer[3, -4, 5, 4])) == Integer[3, -4, 5]
    @test @inferred(unique(iseven, Integer[3, -4, 5, 4]; seen=Set{Bool}())) == Integer[3, -4]
    @test @inferred(unique(n -> n % 3, [5, 1, 8, 9, 3, 4, 10, 7, 2, 6])) == [5, 1, 9]
    for r = (Base.OneTo(-1), Base.OneTo(0), Base.OneTo(1), Base.OneTo(5),
             1:0, 1:1, 1:2, 1:10, 1:.5:.5, 1:.5:1, 1:.5:10, 3:-2:5, 3:-2:3, 3:-2:1,
             StepRangeLen(1.0, 2.0, 0), StepRangeLen(1.0, 2.0, 2), StepRangeLen(1.0, 2.0, 3),
             StepRangeLen(1.0, 0.0, 0), StepRangeLen(1.0, -0.0, 1), StepRangeLen(1.0, 0.0, 2),
             LinRange(1, 2, 3), LinRange(1, 1, 0), LinRange(1, 1, 1), LinRange(1, 1, 10))
        @test @inferred(unique(r)) == invoke(unique, Tuple{Any}, r)
    end
end

@testset "issue 20105" begin
    @test @inferred(unique(x for x in 1:1)) == [1]
    @test unique(x for x in Any[1, 1.0])::Vector{Real} == [1]
    @test unique(x for x in Real[1, 1.0])::Vector{Real} == [1]
    @test @inferred(unique(Integer[1, 1, 2]))::Vector{Integer} == [1, 2]
    @test unique(x for x in []) isa Vector{Any}
end

@testset "unique!" begin
    u = []
    @test unique!(u) === u
    u = [1,1,3,2,1]
    @inferred(unique!(u))
    @test u == [1,3,2]
    @test @inferred(unique!([])) == []
    @test @inferred(unique!(Float64[])) == Float64[]
    u = [1,2,2,3,5,5]
    @test unique!(u) === u
    @test u == [1,2,3,5]
    u = [6,5,5,3,3,2,1]
    @test unique!(u) === u
    @test u == [6,5,3,2,1]
    u = OffsetArray([1,2,2,3,5,5], -1)
    @test unique!(u) === u
    @test u == OffsetArray([1,2,3,5], -1)
    u = OffsetArray([5,5,4,4,2,2,0,-1,-1], -1)
    @test unique!(u) === u
    @test u == OffsetArray([5,4,2,0,-1], -1)
    u = OffsetArray(["w","we","w",5,"r",5,5], -1)
    @test unique!(u) === u
    @test u == OffsetArray(["w","we",5,"r"], -1)
    u = [0.0,-0.0,1.0,2]
    @test unique!(u) === u
    @test u == [0.0,-0.0,1.0,2.0]
    u = [1,NaN,NaN,3]
    @test unique!(u) === u
    @test u[1] == 1
    @test isnan(u[2])
    @test u[3] == 3
    u = [5,"w","we","w","r",5,"w"]
    unique!(u)
    @test u == [5,"w","we","r"]
    u = [1,2,5,1,3,2]
    @test unique!(x -> x ^ 2, [1, -1, 3, -3, 5, -5]) == [1, 3, 5]
    @test unique!(n -> n % 3, [5, 1, 8, 9, 3, 4, 10, 7, 2, 6]) == [5, 1, 9]
    @test @inferred(unique!(iseven, [2, 3, 5, 7, 9])) == [2, 3]
    @test @inferred(unique!(x -> x % 2 == 0 ? :even : :odd, [1, 2, 3, 4, 2, 2, 1])) == [1, 2]
    @test @inferred(unique!(x -> x % 2 == 0 ? :even : "odd", [1, 2, 3, 4, 2, 2, 1]; seen=Set{Union{Symbol,String}}())) == [1, 2]

    @test isempty(unique!(Union{}[]))
    @test eltype(unique!([i for i in ["1"] if i isa Int])) <: Union{}
end

@testset "allunique" begin
    @test allunique([])
    @test allunique(Set())
    @test allunique([1,2,3])
    @test allunique([1 2; 3 4])
    @test allunique([:a,:b,:c])
    @test allunique(Set([1,2,3]))
    @test !allunique([1,1,2])
    @test !allunique([:a,:b,:c,:a])
    @test allunique(unique(randn(100)))  # longer than 32
    @test allunique(collect('A':'z')) # 58-element Vector{Char}
    @test !allunique(repeat(1:99, 1, 2))
    @test !allunique(vcat(pi, randn(1998), pi))  # longer than 1000
    @test allunique(eachrow(hcat(1:10, 1:10)))
    @test allunique(x for x in 'A':'Z' if randn()>0)
    @test !allunique(x for x in repeat(1:2000, 3) if true)
    @test allunique([0.0, -0.0])
    @test allunique(x for x in [0.0, -0.0] if true)
    @test !allunique([NaN, NaN])
    @test !allunique(x for x in [NaN, NaN] if true)
    # ranges
    @test allunique(4:7)
    @test allunique(1:1)
    @test allunique(4.0:0.3:7.0)
    @test allunique(4:-1:5)       # empty range
    @test allunique(7:-1:1)       # negative step
    @test allunique(Date(2018, 8, 7):Day(1):Date(2018, 8, 11))  # JuliaCon 2018
    @test allunique(DateTime(2018, 8, 7):Hour(1):DateTime(2018, 8, 11))
    @test allunique(('a':1:'c')[1:2]) == true
    @test allunique(collect(1:1001))
    for r = (Base.OneTo(-1), Base.OneTo(0), Base.OneTo(1), Base.OneTo(5),
             1:0, 1:1, 1:2, 1:10, 1:.5:.5, 1:.5:1, 1:.5:10, 3:-2:5, 3:-2:3, 3:-2:1,
             StepRangeLen(1.0, 2.0, 0), StepRangeLen(1.0, 2.0, 2), StepRangeLen(1.0, 2.0, 3),
             StepRangeLen(1.0, 0.0, 0), StepRangeLen(1.0, -0.0, 1), StepRangeLen(1.0, 0.0, 2),
             LinRange(1, 2, 3), LinRange(1, 1, 0), LinRange(1, 1, 1), LinRange(1, 1, 10))
        @test allunique(r) == invoke(allunique, Tuple{Any}, r)
    end
    # tuples
    @test allunique(())
    @test allunique((1,2,3))
    @test allunique(ntuple(identity, 40))
    @test !allunique((1,2,3,4,3))
    @test allunique((0.0, -0.0))
    @test !allunique((NaN, NaN))
end

@testset "allequal" begin
    @test allequal(Set())
    @test allequal(Set(1))
    @test !allequal(Set([1, 2]))
    @test allequal(Dict())
    @test allequal(Dict(:a => 1))
    @test !allequal(Dict(:a => 1, :b => 2))
    @test allequal([])
    @test allequal([1])
    @test allequal([1, 1])
    @test !allequal([1, 1, 2])
    @test allequal([:a, :a])
    @test !allequal([:a, :b])
    @test !allequal(1:2)
    @test allequal(1:1)
    @test !allequal(4.0:0.3:7.0)
    @test allequal(4:-1:5)       # empty range
    @test !allequal(7:-1:1)       # negative step
    @test !allequal(Date(2018, 8, 7):Day(1):Date(2018, 8, 11))  # JuliaCon 2018
    @test !allequal(DateTime(2018, 8, 7):Hour(1):DateTime(2018, 8, 11))
    @test allequal(StepRangeLen(1.0, 0.0, 2))
    @test !allequal(StepRangeLen(1.0, 1.0, 2))
    @test allequal(LinRange(1, 1, 0))
    @test allequal(LinRange(1, 1, 1))
    @test allequal(LinRange(1, 1, 2))
    @test !allequal(LinRange(1, 2, 2))
end

@testset "filter(f, ::$S)" for S = (Set, BitSet)
    s = S([1,2,3,4])
    @test s !== filter( isodd, s) == S([1,3])
    @test s === filter!(isodd, s) == S([1,3])
end
@testset "first" begin
    @test_throws ArgumentError first(Set())
    @test first(Set(2)) == 2
end
@testset "pop!" begin
    s = Set(1:5)
    @test 2 in s
    @test pop!(s, 2) == 2
    @test !(2 in s)
    @test_throws KeyError pop!(s, 2)
    @test pop!(s, 2, ()) == ()
    @test 3 in s
    @test pop!(s, 3, ()) == 3
    @test !(3 in s)
    @test pop!(Set(1:2), 2, nothing) == 2
end

@testset "convert" begin
    iset = Set([17, 4711])
    cfset = convert(Set{Float64}, iset)
    @test typeof(cfset) == Set{Float64}
    @test cfset == iset
    fset = Set([17.0, 4711.0])
    ciset = convert(Set{Int}, fset)
    @test typeof(ciset) == Set{Int}
    @test ciset == fset
    ssset = Set(split("foo bar"))
    cssset = convert(Set{String}, ssset)
    @test typeof(cssset) == Set{String}
    @test cssset == Set(["foo", "bar"])
end

@testset "fuzzy testing Set & BitSet" begin
    b1, b2 = rand(-1000:1000, 2)
    e1 = rand(b1-9:1000) # -9 to have an empty list sometimes
    e2 = rand(b2-9:1000)
    l1, l2 = rand(1:1000, 2)
    a1 = b1 <= e1 ? rand(b1:e1, l1) : Int[]
    a2 = b2 <= e2 ? rand(b2:e2, l2) : Int[]
    s1, s2 = Set(a1), Set(a2)
    t1, t2 = BitSet(a1), BitSet(a2)

    for (s, t) = ((s1, t1), (s2, t2))
        @test length(s) == length(t)
        @test issubset(s, t)
        @test issubset(t, s)
        @test isempty(s) == isempty(t)
        isempty(s) && continue
        @test maximum(s) == maximum(t)
        @test minimum(s) == minimum(t)
        @test extrema(s) == extrema(t)
        rs, rt = rand(s), rand(t)
        @test rs in s
        @test rt in s
        @test rs in t
        @test rt in t
        for y in (rs, rt)
            ss = copy(s)
            tt = copy(t)
            pop!(ss, y)
            pop!(tt, y)
            @test BitSet(ss) == tt
            @test Set(tt) == ss
            z = rand(1001:1100) # z ∉ s or t
            push!(ss, z)
            push!(tt, z)
            @test BitSet(ss) == tt
            @test Set(tt) == ss
        end
    end

    res = Dict{String,Union{Bool,Vector{Int}}}()
    function check(desc, val)
        n = val isa Bool ? val : sort!(collect(val))
        r = get!(res, desc, n)
        if n isa Bool || r !== n
            @test r == n
        end
    end
    asbitset(x) = x isa BitSet ? x : BitSet(x)
    asset(x) = x isa Set ? x : Set(x)

    for x1 = (s1, t1), x2 = (s2, t2)
        check("union", union(x1, x2))
        check("intersect", intersect(x1, x2))
        check("symdiff", symdiff(x1, x2))
        check("setdiff", setdiff(x1, x2))
        check("== as Bitset", asbitset(x1) == asbitset(x2))
        check("== as Set", asset(x1) == asset(x2))
        check("issubset", issubset(x1, x2))
        if typeof(x1) == typeof(x2)
            check("<", x1 < x2)
            check("<=", x1 > x2)
            check("union!", union!(copy(x1), x2))
            check("setdiff!", setdiff!(copy(x1), x2))
            x1 isa Set && continue
            check("intersect!", intersect!(copy(x1), x2))
            check("symdiff!", symdiff!(copy(x1), x2))
        end
    end
end

@testset "replace! & replace" begin
    a = [1, 2, 3, 1]
    @test replace(x -> iseven(x) ? 2x : x, a) == [1, 4, 3, 1]
    @test replace(x -> iseven(x) ? 2x : x, Tuple(a)) === (1, 4, 3, 1)
    @test replace!(x -> iseven(x) ? 2x : x, a) === a
    @test a == [1, 4, 3, 1]
    @test replace(a, 1=>0) == [0, 4, 3, 0]
    @test replace(Tuple(a), 1=>0) === (0, 4, 3, 0)
    for count = (1, 0x1, big(1))
        @test replace(a, 1=>0, count=count) == [0, 4, 3, 1]
        @test replace(Tuple(a), 1=>0, count=count) === (0, 4, 3, 1)
    end
    @test replace!(a, 1=>2) === a
    @test a == [2, 4, 3, 2]
    @test replace!(x->2x, a, count=0x2) == [4, 8, 3, 2]

    d = Dict(1=>2, 3=>4)
    @test replace!(x -> x.first > 2 ? x.first=>2*x.second : x, d) === d
    @test d == Dict(1=>2, 3=>8)
    @test replace(d, (3=>8)=>(0=>0)) == Dict(1=>2, 0=>0)
    @test replace!(d, (3=>8)=>(2=>2)) === d
    @test d == Dict(1=>2, 2=>2)

    s = Set([1, 2, 3])
    @test replace(x -> x > 1 ? 2x : x, s) == Set([1, 4, 6])
    for count = (1, 0x1, big(1))
        @test replace(x -> x > 1 ? 2x : x, s, count=count) in [Set([1, 4, 3]), Set([1, 2, 6])]
    end
    @test replace(s, 1=>4) == Set([2, 3, 4])
    @test replace!(s, 1=>2) === s
    @test s == Set([2, 3])
    @test replace!(x->2x, s, count=0x1) in [Set([4, 3]), Set([2, 6])]

    for count = (0, 0x0, big(0)) # count == 0 --> no replacements
        @test replace([1, 2], 1=>0, 2=>0; count) == [1, 2]
        @test replace((1, 2), 1=>0, 2=>0; count) === (1, 2)
        for dict = (Dict(1=>2, 2=>3), IdDict(1=>2, 2=>3))
            @test replace(dict, (1=>2) => (1=>3); count) == dict
        end
        @test replace(Set([1, 2]), 2=>-1; count) == Set([1, 2])
    end

    # test collisions with AbstractSet/AbstractDict
    @test replace!(x->2x, Set([3, 6])) == Set([6, 12])
    @test replace!(x->2x, Set([1:20;])) == Set([2:2:40;])
    @test replace!(kv -> (2kv[1] => kv[2]), Dict(1=>2, 2=>4, 4=>8, 8=>16)) == Dict(2=>2, 4=>4, 8=>8, 16=>16)

    # avoid recursive call issue #25384
    @test_throws MethodError replace!("")

    # test eltype promotion
    x = @inferred replace([1, 2], 2=>2.5)
    @test x == [1, 2.5] && x isa Vector{Float64}

    x = @inferred replace([1, 2], 2=>missing)
    @test isequal(x, [1, missing]) && x isa Vector{Union{Int, Missing}}

    @test_broken @inferred replace([1, missing], missing=>2)
    x = replace([1, missing], missing=>2)
    @test x == [1, 2] && x isa Vector{Int}
    x = @inferred replace([1, missing], missing=>2, count=1)
    @test x == [1, 2] && x isa Vector{Union{Int, Missing}}
    x = @inferred replace([1, missing], missing=>missing)
    @test isequal(x, [1, missing]) && x isa Vector{Union{Int, Missing}}
    x = @inferred replace([1, missing], missing=>2, 1=>missing)
    @test isequal(x, [missing, 2]) && x isa Vector{Union{Int, Missing}}

    # eltype promotion for dicts
    d = Dict(1=>2, 3=>4)
    f = replace(d, (1=>2) => (1=>nothing))
    @test f == Dict(3=>4, 1=>nothing)
    @test eltype(f) == Pair{Int, Union{Nothing, Int}}
    f = replace(d, (1=>2) => (1=>missing), (3=>4)=>(3=>missing))
    @test valtype(f) == Union{Missing,Int}
    f = replace(d, (1=>2) => (1=>'a'), (3=>4)=>(3=>'b'))
    @test valtype(f) == Any
    @test f == Dict(3=>'b', 1=>'a')

    # eltype promotion for sets
    s = Set([1, 2, 3])
    f = replace(s, 2=>missing, 3=>nothing)
    @test f == Set([1, missing, nothing])
    @test eltype(f) == Union{Int,Missing,Nothing}
    f = replace(s, 2=>'a')
    @test eltype(f) == Any
    @test f == Set([1, 3, 'a'])

    # test that isequal is used
    @test replace([NaN, 1.0], NaN=>0.0) == [0.0, 1.0]
    @test replace((NaN, 1.0), NaN=>0.0) === (0.0, 1.0)
    @test replace([1, missing], missing=>0) == [1, 0]
    @test replace((1, missing), missing=>0) === (1, 0)
end

@testset "⊆, ⊊, ⊈, ⊇, ⊋, ⊉, <, <=, issetequal" begin
    a = [2, 1, 2]
    b = [2, 3, 1, 3]
    ua = unique(a)
    ub = unique(b)
    for TA in (Tuple, identity, Set, BitSet, Base.IdSet{Int}),
        TB in (Tuple, identity, Set, BitSet, Base.IdSet{Int}),
        uA = false:true,
        uB = false:true
        A = TA(uA ? ua : a)
        B = TB(uB ? ub : b)
        @test A ⊆ B
        @test A ⊊ B
        @test !(A ⊈ B)
        @test !(A ⊇ B)
        @test !(A ⊋ B)
        @test A ⊉ B
        @test !(B ⊆ A)
        @test !(B ⊊ A)
        @test B ⊈ A
        @test B ⊇ A
        @test B ⊋ A
        @test !(B ⊉ A)
        @test !issetequal(A, B)
        @test !issetequal(B, A)
        for T = (Tuple, identity, Set, BitSet, Base.IdSet{Int})
            @test issetequal(A, T(A))
            @test issetequal(B, T(B))
        end
        if A isa AbstractSet && B isa AbstractSet
            @test A <= B
            @test A <  B
            @test !(A >= B)
            @test !(A >  B)
            @test !(B <= A)
            @test !(B <  A)
            @test B >= A
            @test B >  A
        end
    end
    # first doesn't have length
    @test issetequal(Iterators.filter(iseven, 1:10), [2, 4, 6, 8, 10])
    # both don't have length
    @test issetequal(Iterators.filter(iseven, 1:10), Iterators.filter(iseven, 1:10))
end

@testset "optimized union! with max_values" begin
    # issue #30315
    T = Union{Nothing, Bool}
    @test Base.max_values(T) == 3
    d = Set{T}()
    union!(d, (nothing, true, false))
    @test length(d) == 3
    @test d == Set((nothing, true, false))
    @test nothing in d
    @test true    in d
    @test false   in d

    for X = (Int8, Int16, Int32, Int64)
        @test Base.max_values(Union{Nothing, X}) == (sizeof(X) < sizeof(Int) ?
                                                     2^(8*sizeof(X)) + 1 :
                                                     typemax(Int))
    end
    # this does not account for non-empty intersections of the unioned types
    @test Base.max_values(Union{Int8,Int16}) == 2^8 + 2^16
end

struct OpenInterval{T}
    lower::T
    upper::T
end
Base.in(x, i::OpenInterval) = i.lower < x < i.upper
Base.IteratorSize(::Type{<:OpenInterval}) = Base.SizeUnknown()

@testset "Continuous sets" begin
    i = OpenInterval(2, 4)
    @test 3 ∈ i
    @test issubset(3, i)
end

@testset "IdSet" begin
    a = [1]
    b = [2]
    c = [3]
    d = [4]
    e = [5]
    A = Base.IdSet{Vector{Int}}([a, b, c, d])
    @test !isempty(A)
    B = copy(A)
    @test A ⊆ B
    @test B ⊆ A
    A = filter!(x->isodd(x[1]), A)
    @test A ⊆ B
    @test !(B ⊆ A)
    @test !isempty(A)
    a_ = pop!(A, a)
    @test a_ === a
    @test !isempty(A)
    e_ = pop!(A, a, e)
    @test e_ === e
    @test !isempty(A)
    A = empty!(A)
    @test isempty(A)
end

@testset "⊊, ⊋" begin
    @test !((1, 2) ⊊ (1, 2, 2))
    @test !((1, 2, 2) ⊋ (1, 2))
end

@testset "AbstractSet & Fallback" begin
    mutable struct TestSet{T} <: AbstractSet{T}
        set::Set{T}
        function TestSet{T}() where T
            new{T}(Set{T}())
        end
    end
    set = TestSet{Any}()
    @test sizehint!(set, 1) === set
end

max_stars_repo_name: Physolia/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Nexuscompute/julia


text:
module TestStatistics
    using NullableStats
    using Base.Test

    srand(1)

    for N in (10, 100)
        A = rand(N)
        M = rand(Bool, N)
        mu_A = mean(A)
        nmu_A = Nullable(mu_A)
        i = rand(1:N)
        M[i] = true
        j = rand(1:N)
        while j == i
            j = rand(1:N)
        end
        M[j] = false
        X = NullableArray(A)
        Y = NullableArray(A, M)
        J = find(x -> !x, M)
        B = A[J]
        mu_B = mean(B)
        nmu_B = Nullable(mu_B)

        C = rand(N)
        V = WeightVec(C)
        R = rand(Bool, N)
        R[rand(1:N)] = true
        R[j] = false
        # W = WeightVec(NullableArray(C, R))

        K = find(x -> !x, R)
        L = find(i -> (!M[i] & !R[i]), [1:N...])

        # For testing Base.varzm
        D1 = rand(round(Int, N / 2))
        D2 = -1 .* D1
        D = [D1; D2]
        while mean(D) != 0
            D1 = rand(round(Int, N / 2))
            D2 = -1 .* D1
            D = [D1; D2]
        end
        Q = NullableArray(D)
        S = rand(Bool, round(Int, N / 2))
        U = NullableArray(D, [S; S])
        E = [D[find(x->!x, S)]; D[find(x->!x, S)]]

        @test_throws NullException varm(Y, Nullable{Float64}())

        # Test mean
        for skip in (true, false)
            v = mean(X; skipnull=skip)
            @test_approx_eq v.value mean(A)
            @test !v.isnull

            v = mean(Y; skipnull=skip)
            if skip == false
                @test isequal(v, Nullable{Float64}())
            else
                @test_approx_eq v.value mean(B)
                @test !v.isnull
            end

            v = mean(X, V; skipnull=skip)
            @test_approx_eq v.value mean(A, V)
            @test !v.isnull

            v = mean(Y, V; skipnull=skip)
            if skip == false
                @test isequal(v, Nullable{Float64}())
            else
                @test_approx_eq v.value mean(B, WeightVec(C[J]))
                @test !v.isnull
            end

            # Following tests need to wait until WeightVec constructor is
            # implemented for NullableArray argument
            #
            # v = mean(X, W, skipnull=skip)
            # if skip == false
            #     @test isequal(v, Nullable{Float64}())
            # else
            #     @test_approx_eq v.value mean(A[K], WeightVec(C[K]))
            #     @test !v.isnull
            # end
            # v = mean(Y, W, skipnull=skip)
            # if skip == false
            #     @test isequal(v, Nullable{Float64}())
            # else
            #     @test_approx_eq v.value mean(A[L], WeightVec(C[L]))
            #     @test !v.isnull
            # end
        end

        for corr in (true, false), skip in (true, false)
            # Test Base.varzm
            v = Base.varzm(Q, corrected=corr, skipnull=skip)
            @test_approx_eq v.value Base.varzm(D, corrected=corr)
            @test !v.isnull

            v = Base.varzm(U, corrected=corr, skipnull=skip)
            if skip == false
                @test isequal(v, Nullable{Float64}())
            else
                @test_approx_eq v.value Base.varzm(E, corrected=corr)
                @test !v.isnull
            end

            # Test varm, stdm
            for method in (varm, stdm)
                for mu in (mu_A, nmu_A)
                    v = method(X, mu, corrected=corr, skipnull=skip)
                    @test_approx_eq v.value method(A, mu_A, corrected=corr)
                    @test !v.isnull
                end

                for mu in (mu_B, nmu_B)
                    v = method(Y, mu, corrected=corr, skipnull=skip)
                    if skip == false
                        @test isequal(v, Nullable{Float64}())
                    else
                        @test_approx_eq v.value method(B, mu_B, corrected=corr)
                        @test !v.isnull
                    end
                end
            end

            # Test var, std
            for method in (var, std)
                for mu in (nothing, mu_A, nmu_A, :a)
                    if mu == :a
                        @test_throws ErrorException method(X, mean=mu,
                                                               corrected=corr,
                                                               skipnull=skip)
                    else
                        v = method(X, mean=mu, corrected=corr, skipnull=skip)
                        @test_approx_eq v.value method(A, mean=mu_A, corrected=corr)
                        @test !v.isnull
                    end
                end

                for mu in (0, Nullable(0))
                    v = method(X, mean=mu, corrected=corr, skipnull=skip)
                    @test_approx_eq v.value method(A, mean=0, corrected=corr)
                    @test !v.isnull
                end

                for mu in (nothing, mu_B, nmu_B)
                    v = method(Y, mean=mu, corrected=corr, skipnull=skip)
                    if skip == false
                        @test isequal(v, Nullable{Float64}())
                    else
                        @test_approx_eq v.value method(B, mean=mu_B, corrected=corr)
                        @test !v.isnull
                    end
                end

                for mu in (0, Nullable(0))
                    v = method(Y, mean=mu, corrected=corr, skipnull=skip)
                    if skip == false
                        @test isequal(v, Nullable{Float64}())
                    else
                        @test_approx_eq v.value method(B, mean=0, corrected=corr)
                        @test !v.isnull
                    end
                end
            end
        end
    end
end

max_stars_repo_name: davidagold/NullableStats

max_issues_repo_name: davidagold/NullableStats

max_forks_repo_name: davidagold/NullableStats


text:
"""
Memory-Mapped implementation of Kernel Spectral Clustering
"""
function kscbignet(network_file::String, file_dir::String, file_name::String, delimiter::Char, maxMB::Int; fraction=0.15::Float64, method="furs"::String, decomp="eigs"::String, minK=2::Int, maxK=100::Int, convertF=true::Bool)
  if minK < 2 || minK > maxK
    error("minK must be at least 2 and smaller than maxK")
  end
  if !isdir("$(file_dir)/$(file_name)")
    mkdir("$(file_dir)/$(file_name)")
  end
  resultfile = open("$(file_dir)/$(file_name)/resultKSC.txt", "w")
  if convertF
    println("File Conversion $(network_file)")
    tic()
    fileCount, firstNode, lastNode, fileWeighted = convertDSVNet(network_file, file_dir, file_name, delimiter, maxMB)
    timeConversion = toq()
    println("File Conversion step finished, time elapsed $(timeConversion)s")
    write(resultfile, "File Conversion step: time $(timeConversion)s\n")
  else
    file = open("$(file_dir)/$(file_name)/info.bin")
    fileCount = read(file, Int)
    firstNode = read(file, Int)
    lastNode = read(file, Int)
    fileWeighted = read(file, Int)
    close(file)
  end
  println("Assigning global variables")
  # stores mmaps of nodeArray, degreeArray, neighbourMatrix, weightMatrix for each file on each processor
  for p in procs()
    if fileWeighted == 0
      @spawnat(p, netmapU(file_dir, file_name, fileCount))
    else
      @spawnat(p, netmapW(file_dir, file_name, fileCount))
    end
    @spawnat(p, qmap("$(file_dir)/$(file_name)/qtest.bin", firstNode, lastNode))
  end
  # TODO make max subset size dynamic based on maxMB?
  sampleSize = min(5000, int(ceil((lastNode-firstNode)*fraction)))
  sampleSizeBlock = int(floor(sampleSize/fileCount))
  println("Total sample size for $(fileCount) files = $(sampleSize)")
  println("Subset selection")
  tic()
  # perform the subset selection on each file in parallel (map) and concatenate them together (reduce)
  if fileWeighted == 0
    trainNodes, trainNeighbours, validNodes, validNeighbours = @parallel (subsetCombineU) for fileIndex = 1:fileCount
      # subset selection degree threshold
      degreeThreshold = mean(fileData[fileIndex][2])
      trainIndices, trainComplete, validIndices, validComplete = selectSubset(method, fileIndex, sampleSizeBlock, sampleSizeBlock, degreeThreshold)
      if !trainComplete || !validComplete
        println("Incomplete subset selection on file set $(fileIndex), retrying with median degree threshold")
        degreeThreshold = median(fileData[fileIndex][2])
        trainIndices, trainComplete, validIndices, validComplete = selectSubset(method, fileIndex, sampleSizeBlock, sampleSizeBlock, degreeThreshold)
        if !trainComplete || !validComplete
          println("Incomplete subset selection on file set $(fileIndex), continuing")
        end
      end
      trainNodeSeq = Array(Int, 0)
      trainNeighbourSeq = Array(Int, 0)
      validNodeSeq = Array(Int, 0)
      validNeighbourSeq = Array(Int, 0)
      generateSubsetU(fileIndex, trainIndices, trainNodeSeq, trainNeighbourSeq)
      generateSubsetU(fileIndex, validIndices, validNodeSeq, validNeighbourSeq)
      # return training and validation data as tuple
      (trainNodeSeq, trainNeighbourSeq, validNodeSeq, validNeighbourSeq)
    end
  else
    trainNodes, trainNeighbours, trainWeights, validNodes, validNeighbours, validWeights = @parallel (subsetCombineW) for fileIndex = 1:fileCount
      # subset selection degree threshold
      degreeThreshold = mean(fileData[fileIndex][2])
      trainIndices, trainComplete, validIndices, validComplete = selectSubset(method, fileIndex, sampleSizeBlock, sampleSizeBlock, degreeThreshold)
      if !trainComplete || !validComplete
        println("Incomplete subset selection on file set $(fileIndex), retrying with median degree threshold")
        degreeThreshold = median(fileData[fileIndex][2])
        trainIndices, trainComplete, validIndices, validComplete = selectSubset(method, fileIndex, sampleSizeBlock, sampleSizeBlock, degreeThreshold)
        if !trainComplete || !validComplete
          println("Incomplete subset selection on file set $(fileIndex), continuing")
        end
      end
      trainNodeSeq = Array(Int, 0)
      trainNeighbourSeq = Array(Int, 0)
      trainWeightSeq = Array(Float64, 0)
      validNodeSeq = Array(Int, 0)
      validNeighbourSeq = Array(Int, 0)
      validWeightSeq = Array(Float64, 0)
      generateSubsetW(fileIndex, trainIndices, trainNodeSeq, trainNeighbourSeq, trainWeightSeq)
      generateSubsetW(fileIndex, validIndices, validNodeSeq, validNeighbourSeq, validWeightSeq)
      # return training and validation data as tuple
      (trainNodeSeq, trainNeighbourSeq, trainWeightSeq, validNodeSeq, validNeighbourSeq, validWeightSeq)
    end
  end
  timeSubset = toq()
  println("Subset selection step finished, time elapsed $(timeSubset)s")
  # TODO: store to file?
  trainSubset = unique(trainNodes)
  validSubset = unique(validNodes)
  write(resultfile, "Subset selection step: sizeT $(length(trainSubset)), sizeV $(length(validSubset)), time $(timeSubset)s\n")
  # change node indices to be sequential to avoid zeros in matrix representation
  permIndices = trainSubset
  trainN = length(permIndices)
  #[trainNodes[trainNodes .== permIndices[value]] = value for value in 1:length(permIndices)]
  @inbounds for i in 1:length(trainNodes)
    trainNodes[i] = findfirst(permIndices, trainNodes[i])
  end
  permIndices = validSubset
  validN = length(permIndices)
  #[validNodes[validNodes .== permIndices[value]] = value for value in 1:length(permIndices)]
  @inbounds for i in 1:length(validNodes)
    validNodes[i] = findfirst(permIndices, validNodes[i])
  end
  permIndices = []
  # avoid zero index in neighbours
  lastNeighbour = lastNode
  if firstNode == 0
    lastNeighbour += 1
    trainNeighbours[trainNeighbours .== 0] = lastNeighbour
    validNeighbours[validNeighbours .== 0] = lastNeighbour
  end
  # training phase
  println("Training model")
  tic()
  if fileWeighted == 0
    modelSparse = sparse(trainNodes, trainNeighbours, fill(1.0, length(trainNodes)), trainN, lastNeighbour)
  else
    modelSparse = sparse(trainNodes, trainNeighbours, trainWeights, trainN, lastNeighbour)
  end
  trainNorm, λ, α, β, CB = modelTrain(modelSparse, trainN, maxK, decomp)
  timeTrain = toq()
  println("Training phase finished, time elapsed $(timeTrain)s")
  write(resultfile, "Training phase: time $(timeTrain)s\n")
  # validation phase
  if fileWeighted == 0
    modelSparse = sparse(validNodes, validNeighbours, fill(1.0, length(validNodes)), validN, lastNeighbour)
  else
    modelSparse = sparse(validNodes, validNeighbours, validWeights, validN, lastNeighbour)
  end
  println("Validating model")
  tic()
  baf, kBAF, indexBAF, valueBAF = modelValid(modelSparse, validN, trainNorm, maxK, minK, α, β, CB)
  timeValid = toq()
  println("Validation phase finished (k = $(kBAF)), time elapsed $(timeValid)s")
  write(resultfile, "Validation phase: BAF $(valueBAF), k $(kBAF), time $(timeValid)s\n")
  modelSparse = []
  # test phase
  println("Test phase")
  tic()
  if fileWeighted == 0
    modelTestU(fileCount, trainNorm, firstNode, lastNode, α, β, CB, kBAF, maxMB)
  else
    modelTestW(fileCount, trainNorm, firstNode, lastNode, α, β, CB, kBAF, maxMB)
  end
  timeTest = toq()
  println("Test phase finished, time elapsed $(timeTest)s")
  write(resultfile, "Test phase: time $(timeTest)s\n")
  close(resultfile)
  return qData, baf, kBAF, firstNode, lastNode, trainSubset, fileCount, fileWeighted
end

function subsetCombineU(x::Tuple, y::Tuple)
  return vcat(x[1], y[1]), vcat(x[2], y[2]), vcat(x[3], y[3]), vcat(x[4], y[4])
end

function subsetCombineW(x::Tuple, y::Tuple)
  return vcat(x[1], y[1]), vcat(x[2], y[2]), vcat(x[3], y[3]), vcat(x[4], y[4]), vcat(x[5], y[5]), vcat(x[6], y[6])
end

function generateSubsetU(fileIndex::Int, indicesS::Array{Int}, nodesS::Array{Int}, neighboursS::Array{Int})
  for i in indicesS
    tempNeighbourSeq = fileData[fileIndex][3][:, i]
    # remove the -1 and 0 padding
    foundIndex = findfirst(tempNeighbourSeq, -1)
    if foundIndex != 0
      lengthSeq = length(tempNeighbourSeq)
      splice!(tempNeighbourSeq, foundIndex:lengthSeq)
    end
    # remove links between training nodes
    for j in indicesS
      foundIndex = findfirst(tempNeighbourSeq, j)
      if foundIndex != 0
        splice!(tempNeighbourSeq, foundIndex)
      end
    end
    push!(nodesS, fill(fileData[fileIndex][1][i], length(tempNeighbourSeq))...)
    push!(neighboursS, tempNeighbourSeq...)
  end
end

function generateSubsetW(fileIndex::Int, indicesS::Array{Int}, nodesS::Array{Int}, neighboursS::Array{Int}, weightsS::Array{Float64})
  for i in indicesS
    tempNeighbourSeq = fileData[fileIndex][3][:, i]
    tempWeightSeq = fileData[fileIndex][4][:, i]
    # remove the -1 and 0 padding
    foundIndex = findfirst(tempNeighbourSeq, -1)
    if foundIndex != 0
      lengthSeq = length(tempNeighbourSeq)
      splice!(tempNeighbourSeq, foundIndex:lengthSeq)
      splice!(tempWeightSeq, foundIndex:lengthSeq)
    end
    # remove links between training nodes
    for j in indicesS
      foundIndex = findfirst(tempNeighbourSeq, j)
      if foundIndex != 0
        splice!(tempNeighbourSeq, foundIndex)
        splice!(tempWeightSeq, foundIndex)
      end
    end
    push!(nodesS, fill(fileData[fileIndex][1][i], length(tempNeighbourSeq))...)
    push!(neighboursS, tempNeighbourSeq...)
    push!(weightsS, tempWeightSeq...)
  end
end

function netmapU(file_dir::String, file_name, fileCount::Int)
  global fileData = Array((Array{Int, 2}, Array{Float64, 2}, Array{Int, 2}), fileCount)
  @inbounds for i = 1:fileCount
    fileData[i] = mmapNetU(file_dir, file_name, i)
  end
end

function netmapW(file_dir::String, file_name, fileCount::Int)
  global fileData = Array((Array{Int, 2}, Array{Float64, 2}, Array{Int, 2}, Array{Float64, 2}), fileCount)
  @inbounds for i = 1:fileCount
    fileData[i] = mmapNetW(file_dir, file_name, i)
  end
end

function qmap(qFile::String, minN::Int, maxN::Int)
  global qData
  file = open(qFile, "w+")
  qData = mmap_array(Int, (maxN-minN+1, 2), file)
  close(file)
end


text:
module WeightedSumPlatforms
using ..Platforms, ..ParameterPaths
import ..Platforms: SolverStyle, dictionary, measure, dualdictionary, param_first,
    param_double, correctparamformat, unsafe_dictionary, param_inbetween
using BasisFunctions: Dictionary, AbstractMeasure, MultiDict

export WeightedSumPlatform
"""
A WeightedSumPlatform is the union of a finite number of copies of a single frame,
each weighted by a function.
"""
struct WeightedSumPlatform <: FramePlatform
    P       :: Platform
    weights :: NTuple{N,Function} where N
end

correctparamformat(platform::WeightedSumPlatform, param::NTuple) =
    all(map(parami->correctparamformat(platform.P, parami), param))

unsafe_dictionary(platform::WeightedSumPlatform, i::NTuple) =
    MultiDict([weight(platform, j) * unsafe_dictionary(platform.P, i[j]) for j in 1:length(platform.weights)])

WeightedSumPlatform(platform::Platform, weights::Function...) where N = WeightedSumPlatform(platform, weights)

SolverStyle(platform::WeightedSumPlatform, ::SamplingStyle) = AZStyle()

weight(platform::WeightedSumPlatform, i) = platform.weights[i]

# dictionary(platform::WeightedSumPlatform, i) =
#     MultiDict([weight(platform, j) * dictionary(platform.P, i) for j in 1:length(platform.weights)])

function dualdictionary(platform::WeightedSumPlatform, param::NTuple, measure::AbstractMeasure; options...)
    denom = (x...)->sum(map(w->abs(w(x...))^2, platform.weights))
    # TODO: discuss, what is the relation between param, L of a platform and platform.P
    MultiDict([((x...)->(platform.weights[j](x...)/denom(x...))) * dualdictionary(platform.P, param[j], measure; options...) for j=1:length(platform.weights)])
end

measure(platform::WeightedSumPlatform) = measure(platform.P)

param_first(platform::WeightedSumPlatform) = ntuple(k->param_first(platform.P),Val(length(platform.weights)))
param_double(platform::WeightedSumPlatform, param::NTuple) =
    ntuple(k->param_double(platform.P, param[k]), Val(length(platform.weights)))
param_inbetween(platform::WeightedSumPlatform, param1::NTuple, param2::NTuple) =
    ntuple(k->param_inbetween(platform.P, param1[k], param2[k]), Val(length(platform.weights)))

using ..ParameterPaths:  default_param_path, NTupleParameterPath, HierarchyPath
import ..ParameterPaths: default_param_path

function default_param_path(platform::WeightedSumPlatform)
    firstpath = default_param_path(platform.P)
    secondpath = NTupleParameterPath{length(platform.weights)}(first(firstpath))
    HierarchyPath(firstpath, secondpath)
end

end


text:
__precompile__()

module LSODA

using Compat, DiffEqBase
import DiffEqBase: solve
using LinearAlgebra, Printf
using LSODA_jll

const warnkeywords =
    (:save_idxs, :d_discontinuities, :isoutofdomain, :unstable_check,
     :calck, :progress, :timeseries_steps, :dense,
     :dtmin, :dtmax,
     :internalnorm, :gamma, :beta1, :beta2, :qmax, :qmin, :qoldinit)

function __init__()
    global warnlist = Set(warnkeywords)
end

abstract type LSODAAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
abstract type LSODAAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
struct lsoda <: LSODAAlgorithm end

export lsoda, lsoda_0, lsoda_opt_t, lsoda_context_t, lsoda_prepare, lsoda_reset, lsoda_opt_t, lsoda_free, lsoda_evolve!, UserFunctionAndData

export LSODAAlgorithm, solve

include("types_and_consts.jl")
include("handle.jl")
include("solver.jl")
include("common.jl")

end # module


text:
"""
ULMFiT - Text Classifier
"""

cd(@__DIR__)
include("custom_layers.jl")     # including custome layers created for ULMFiT model
include("utils.jl")     # including utilities
include("fine_tune_lm.jl")        # including useful functions from fine-tuning file
# include("data_loaders.jl")      # including helper functions for IMDB data loading

mutable struct TextClassifier
    vocab::Vector
    rnn_layers::Flux.Chain
    linear_layers::Flux.Chain
end

function TextClassifier(lm::LanguageModel=load_model!(), clsfr_out_sz::Integer=1, clsfr_hidden_sz::Integer=50, clsfr_hidden_drop::Float64=0.0)
    return TextClassifier(
        lm.vocab,
        lm.layers[1:8],
        Chain(
            gpu(PooledDense(length(lm.lstm_layer3.layer.cell.h), clsfr_hidden_sz, relu)),
            gpu(BatchNorm(clsfr_hidden_sz, relu)),
            Dropout(clsfr_hidden_drop),
            gpu(Dense(clsfr_hidden_sz, clsfr_out_sz, sigmoid)),
        )
    )
end

Flux.@treelike TextClassifier

# Forward step for classifier
function forward(classifier::TextClassifier, gen; now_per_pass::Integer=32)
    X = take!(gen)
    l = length(X)
    maxpools, meanpools = [], []
    # Truncated Backprop through time
    for i=1:floor(l/now_per_pass)   # Tracking is swiched off inside this loop
        (i == 1 && l%now_per_pass != 0) ? (last_idx = l%now_per_pass) : (last_idx = now_per_pass)
        H = broadcast(x -> gpu(classifier.rnn_layers[1](indices(x, classifier.vocab, "_unk_"))), X[1:last_idx])
        H = broadcast(x -> Tracker.data(classifier.rnn_layers[2:end](x)), H)
        H = cat(H..., dims=3)
        push!(maxpools, maximum(H, dims=3)[:, :, 1])
        push!(meanpools, (sum(H, dims=3)/size(H, 3))[:, :, 1])
        X = X[last_idx+1:end]
        Flux.truncate!(classifier.rnn_layers)
        GC.gc()
    end
    # last part of the sequecnes in X - Tracking is swiched on
    H = broadcast(x -> classifier.rnn_layers[1](indices(x, classifier.vocab, "_unk_")), X)
    H = classifier.rnn_layers[2:end].(H)
    H = H[end:-1:1]
    append!(H, maxpools)
    append!(H, meanpools)
    H = classifier.linear_layers(H[end:-1:1])
    return H
end

# Loss function for classifier
function loss(classifier::TextClassifier, gen::Channel; now_per_pass::Integer=32)
    H = forward(classifier, gen, now_per_pass=now_per_pass)
    Y = gpu(take!(gen))
    l = crossentropy(H, Y)
    Flux.reset!(classifier.rnn_layers)
    return l
end

# Training of classifier
function train_classifier!(classifier::TextClassifier=TextClassifier(), classes::Integer=1, hidden_layer_size::Integer=50, data_loader::Channel=imdb_classifier_data;
    batchsize::Integer=64, bptt::Integer=70, gradient_clip::Float64=0.25, stlr_cut_frac::Float64=0.1, stlr_ratio::Number=32, stlr_η_max::Float64=0.01, epochs::Integer=1, checkpoint_itvl=5000)
    trainable = []
    append!(trainable, [classifier.rnn_layers[[1, 3, 5, 7]]...])
    append!(trainable, [classifier.linear_layers[[1, 4]]...])
    opts = [ADAM(0.001, (0.7, 0.99)) for i=1:length(trainable)]
    cut = num_of_iters * epochs * stlr_cut_frac
    gpu!.(classifier.rnn_layers)

    for epoch=1:epochs
        gen = data_loader()
        num_of_iters = take!(gen)
        T = num_of_iters-Int(floor((num_of_iters*2)/100))
        set_trigger!.(T, classifier.rnn_layers)
        for iter=1:num_of_iters
            l = loss(classifier, gen, now_per_pass = now_per_pass)

            # Slanted triangular learning rates
            t = iter + (epoch-1)*num_of_iters
            p_frac = (iter < cut) ? iter/cut : (1 - ((iter-cut)/(cut*(1/stlr_cut_frac-1))))
            ηL = stlr_η_max*((1+p_frac*(stlr_ratio-1))/stlr_ratio)

            # Gradual-unfreezing Step with discriminative fine-tuning
            unfreezed_layers, cur_opts = (epoch < length(trainable)) ? (trainable[end-epoch+1:end], opts[end-epoch+1:end]) : (trainable, opts)
            discriminative_step!(unfreezed_layers, ηL, l, gradient_clip, cur_opts)

            # ASGD Step
            asgd_step!.(iter, classifier.rnn_layers)

            reset_masks!.(classifier.rnn_layers)

            println("loss:", l, "epoch: ", epoch, "iteration: ", iter)
        end
    end
end


text:
"""
    Interval(first::String, last::String)

A type to represent ISO8601 intervals.
"""
struct Interval
    first::String
    last::String
end

JSON.lower(i::Interval) = i.first * '/' * i.last


text:
## Specialized solver for sparse symmetric integer linear systems with exact rational solution

include("Modulos.jl")
using .Modulos
import Base.GMP: MPZ

import LinearAlgebra: BlasInt, checknonsingular, LU, tril!, triu!, ipiv2perm,
                      lu, lu!, ldiv!, Factorization, issuccess
using SparseArrays
import SparseArrays: getcolptr

using BigRationals

function rational_lu!(B::SparseMatrixCSC, col_offset, check=true)
    Tf = eltype(B)
    m, n = size(B)
    minmn = min(m, n)
    info = 0
    @inbounds begin
        for k in 1:minmn
            ipiv = getcolptr(B)[k] + col_offset[k]
            piv = nonzeros(B)[ipiv]
            if iszero(piv)
                check && checknonsingular(k-1, Val(false)) # TODO update with Pivot
                return LU{Tf,SparseMatrixCSC{Tf,Int}}(B, collect(1:minmn), convert(BlasInt, k-1))
            end
            Bkkinv = inv(piv)
            @simd for i in ipiv+1:getcolptr(B)[k+1]-1
                nonzeros(B)[i] *= Bkkinv
            end
            for j in k+1:n
                r1 = getcolptr(B)[j]
                r2 = getcolptr(B)[j+1]-1
                r = searchsortedfirst(rowvals(B), k, r1, r2, Base.Forward)
                ((r > r2) || (rowvals(B)[r] != k)) && continue
                Bkj = nonzeros(B)[r]
                for i in ipiv+1:getcolptr(B)[k+1]-1
                    Bik = nonzeros(B)[i]
                    l = i - ipiv

                    while rowvals(B)[l+r] < rowvals(B)[i]
                        r += 1
                    end
                    nonzeros(B)[l+r] -= Bik * Bkj
                end
            end
        end
    end
    check && checknonsingular(info, Val(false))
    return LU{Tf,SparseMatrixCSC{Tf,Int}}(B, Vector{BlasInt}(1:minmn), convert(BlasInt, info))
end

# function lu!(B::SparseMatrixCSC{<:Rational}, ::Val{Pivot} = Val(false);
#                            col_offset, check::Bool = true) where Pivot
function rational_lu!(B::SparseMatrixCSC{BigRational}, col_offset, check::Bool=true)
    Tf = Rational{BigInt}
    m, n = size(B)
    minmn = min(m, n)
    info = 0
    Bkkinv = BigRational()
    tmp = BigRational()
    @inbounds begin
        for k in 1:minmn
            ipiv = getcolptr(B)[k] + col_offset[k]
            piv = nonzeros(B)[ipiv]
            if iszero(piv)
                check && checknonsingular(k-1, Val(false)) # TODO update with Pivot
                return LU{Tf,SparseMatrixCSC{Tf,Int}}(B, collect(1:minmn), convert(BlasInt, k-1))
            end
            BigRationals.MPQ.inv!(Bkkinv, piv)
            @simd for i in ipiv+1:getcolptr(B)[k+1]-1
                BigRationals.MPQ.mul!(nonzeros(B)[i], Bkkinv)
            end
            for j in k+1:n
                r1 = getcolptr(B)[j]
                r2 = getcolptr(B)[j+1]-1
                r = searchsortedfirst(rowvals(B), k, r1, r2, Base.Forward)
                ((r > r2) || (rowvals(B)[r] != k)) && continue
                Bkj = nonzeros(B)[r]
                for i in ipiv+1:getcolptr(B)[k+1]-1
                    Bik = nonzeros(B)[i]
                    l = i - ipiv

                    while rowvals(B)[l+r] < rowvals(B)[i]
                        r += 1
                    end
                    # Base.GMP.MPZ.mul!(tmp, Bik, Bkj)
                    # Base.GMP.MPZ.sub!(nonzeros(B)[l+r], tmp)
                    BigRationals.MPQ.mul!(tmp, Bik, Bkj)
                    BigRationals.MPQ.sub!(nonzeros(B)[l+r], tmp)
                end
            end
        end
    end
    check && checknonsingular(info, Val(false))
    return LU{Tf,SparseMatrixCSC{Tf,Int}}(Tf.(B), Vector{BlasInt}(1:minmn), convert(BlasInt, info))
end

# function lu(A::SparseMatrixCSC{<:Rational}, pivot::Union{Val{false}, Val{true}} = Val(false); check::Bool = true)
function rational_lu(A::SparseMatrixCSC, check::Bool=true, ::Type{Ti}=BigRational) where {Ti}
    Tf = Ti == BigRational ? Rational{BigInt} : Ti

    Base.require_one_based_indexing(A)
    _I, _J, _V = findnz(A)
    I, J, V = issorted(_J) ? (_I, _J, _V) : begin
        _indices = sortperm(_J)
        @inbounds (_I[_indices], _J[_indices], _V[_indices])
    end
    # @inbounds if !issorted(_J)
    #     indices = sortperm(J)
    #     I = I[indices]; J = J[indices]; V = V[indices]
    # end
    isempty(J) && return LU{Tf,SparseMatrixCSC{Tf,Int}}(A, Int[], convert(BlasInt, 0))
    m, n = size(A)
    minmn = min(m, n)
    if J[1] != 1 || I[1] != 1
        check && checknonsingular(1, Val(false)) # TODO update with Pivot
        # return LU{eltype(A), typeof(A)}(A, collect(1:minmn), convert(BlasInt, 1))
        return LU{Tf,SparseMatrixCSC{Tf,Int}}(A, collect(1:minmn), convert(BlasInt, 1))
    end

    col_offset = zeros(Int, minmn) # for each col, index of the pivot element
    idx_cols = [[I[i] for i in getcolptr(A)[col+1]-1:-1:getcolptr(A)[col]] for col in 1:minmn]
    # For each column, indices of the non-zeros elements
    @inbounds for col in 2:minmn
        sort!(idx_cols[col-1]; rev=true)
        # All idx_cols[x] are sorted by decreasing order for x < col
        # @show idx_cols[col]
        for row_j in idx_cols[col]
            row_j >= col && continue
            col_offset[col] += 1
            # @show idx_cols[row_j]
            for row_i in idx_cols[row_j]
                row_i <= row_j && break # Because the row_i are sorted in decreasing order
                # TODO optimize using the fact that the row_i arrive in decreasing order
                # and the fact that idx_cols[col] is initially sorted in decreasing order
                row_i ∈ idx_cols[col] && continue
                push!(idx_cols[col], row_i)
                push!(J, col)
                push!(I, row_i)
                push!(V, 0)
            end
        end
    end
    B = sparse(I, J, Ti.(V)) # TODO update with Pivot
    # lu!(B, col_offset, check)
    rational_lu!(B, col_offset, check)
end

#=
function lu(A::Hermitian{T, <:SparseMatrixCSC{T}}, pivot::Union{Val{false}, Val{true}} = Val(false); check::Bool = true) where {T<:Rational, Pivot}
    lu(ishermitian(A.data) ? A.data : sparse(A), pivot)
end

function Base.getproperty(F::LU{T,<:SparseMatrixCSC{<:Rational{<:Integer}, <:Integer}}, d::Symbol) where T
    m, n = size(F)
    if d === :L
        L = tril!(getfield(F, :factors)[1:m, 1:min(m,n)])
        for i = 1:min(m,n); L[i,i] = one(T); end
        return L
    elseif d === :U
        return triu!(getfield(F, :factors)[1:min(m,n), 1:n])
    elseif d === :p
        return ipiv2perm(getfield(F, :ipiv), m)
    elseif d === :P
        return Matrix{T}(LinearAlgebra.I, m, m)[:,invperm(F.p)]
    else
        getfield(F, d)
    end
end
=#

function forward_substitution!(L::SparseMatrixCSC, b)
    _, n = size(L)
    _, m = size(b)
    @inbounds for col in 1:n
        k = getcolptr(L)[col]
        rowvals(L)[k] == col || checknonsingular(col, Val(true))
        invnzLk = inv(nonzeros(L)[k])
        x = invnzLk .* b[col,:]
        b[col,:] .= x
        #=@simd=# for i in (k+1):getcolptr(L)[col+1]-1
            nzLi = nonzeros(L)[i]
            rvLi = rowvals(L)[i]
            @simd for j in 1:m
                b[rvLi,j] -= nzLi*x[j]
            end
        end
    end
    nothing
end

function backward_substitution!(U::SparseMatrixCSC, b)
    _, n = size(U)
    _, m = size(b)
    @inbounds for col in n:-1:1
        k = getcolptr(U)[col+1]-1
        rowvals(U)[k] == col || checknonsingular(col, Val(true))
        invnzUk = inv(nonzeros(U)[k])
        x = invnzUk .* b[col,:]
        b[col,:] .= x
        #=@simd=# for i in getcolptr(U)[col]:(k-1)
            nzUi = nonzeros(U)[i]
            rvUi = rowvals(U)[i]
            @simd for j in 1:m
                b[rvUi,j] -= nzUi*x[j]
            end
        end
    end
    nothing
end

function linsolve!(F::LU{<:Any,<:AbstractSparseMatrix}, B::Base.StridedVecOrMat)
    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
    BB = similar(B, TFB, size(B))
    copyto!(BB, B)
    m, n = size(F)
    minmn = min(m,n)
    L = tril!(getfield(F, :factors)[1:m, 1:minmn])
    for i = 1:minmn; L[i,i] = 1; end
    forward_substitution!(L, BB)
    x = triu!(getfield(F, :factors)[1:minmn, 1:n])
    backward_substitution!(x, BB)
    return BB
end

#=
function ldiv!(F::LU{<:Any,<:AbstractSparseMatrix}, B::Base.StridedVecOrMat)
    forward_substitution!(F.L, B)
    backward_substitution!(F.U, B)
    return B
end
=#

function rational_solve(::Val{N}, A, Y) where N
    B = rational_lu(A, false)
    if !issuccess(B)
        error("Singular exception while equilibrating. Is the graph connected?")
    end
    Z = linsolve!(B, Rational{BigInt}.(Y))
    return hcat(zeros(Rational{Int128}, N), Rational{Int128}.(Z)')
    # Rational{Int64} is not enough for tep for instance.
end

# function _inner_dixon_p!(Z::Matrix{Rational{T}}, h, x̄, sqh, tmp) where T
#     for j in eachindex(Z)
#         ua = MPZ.set(h)
#         ub = @inbounds x̄[j]
#         va = zero(T)
#         vb = one(T)
#         k = 0
#         while ub >= sqh
#             k += 1
#             # cpua = deepcopy(ua)
#             # cpub = deepcopy(ub)
#             MPZ.tdiv_qr!(tmp, ua, ua, ub)
#             ua, ub = ub, ua
#             # @assert tmp == cpua ÷ cpub
#             # @assert ua == cpub
#             # @assert ub == cpua - tmp * cpub
#             # cpuc = deepcopy(va)
#             if typemin(Clong) < vb < typemax(Clong)
#                 MPZ.mul_si!(tmp, vb % Clong)
#             else
#                 tmp *= vb
#             end
#             flag = signbit(va)
#             va = abs(va)
#             if va < typemax(Culong)
#                 if flag
#                     MPZ.sub_ui!(tmp, va)
#                 else
#                     MPZ.add_ui!(tmp, va)
#                 end
#                 va, vb = vb, T(tmp)
#             else
#                 va, vb = vb, va + tmp
#             end
#             # @assert vb == cpuc + tmp * va
#         end
#         uc, vc = if T === BigInt
#             Base.divgcd(ub, vb)
#         else
#             ud, vd = Base.divgcd(ub, vb)
#             if T !== BigInt
#                 m = typemin(T)
#                 M = typemax(T)
#                 (m < ud < M && m < vd < M) || return false
#             end
#             (ud % T, vd % T)
#         end
#
#         @inbounds Z[j] = (-2*isodd(k)+1) * Base.checked_den(uc, vc)
#         # @show Z[j]
#         # @assert mod((-1)^isodd(k) * ub, h) == mod(vb * x̄[j], h)
#     end
#     return true
# end

function copyuntil(j, oldZ, ::Type{T}) where T
    Z = similar(oldZ, T)
    for i in eachindex(Z)
        i == j && return Z
        Z[i] = oldZ[i]
    end
    error("Invalid failure of _inner_dixon_p!, please report this error")
    return Z # Does not matter
end

function _inner_dixon_p!(indices, Z::Matrix{Rational{T}}, h, x̄, sqh, tmp) where T
    while !isempty(indices)
        j = pop!(indices)
        ua = MPZ.set(h)
        ub = deepcopy(@inbounds x̄[j])
        va = Int128(0)
        vb = Int128(1)
        k = 0
        while ub >= sqh
            k += 1
            # cpua = deepcopy(ua)
            # cpub = deepcopy(ub)
            MPZ.tdiv_qr!(tmp, ua, ua, ub)
            ua, ub = ub, ua
            # @assert tmp == cpua ÷ cpub
            # @assert ua == cpub
            # @assert ub == cpua - tmp * cpub
            # cpuc = deepcopy(va)
            if typemin(Clong) < vb < typemax(Clong)
                MPZ.mul_si!(tmp, vb % Clong)
            else
                tmp *= vb
            end
            flag = signbit(va)
            va = abs(va)
            if va < typemax(Culong)
                if flag
                    MPZ.sub_ui!(tmp, va)
                else
                    MPZ.add_ui!(tmp, va)
                end
                va, vb = vb, Int128(tmp)
            else
                va, vb = vb, va + tmp
            end
            # @assert vb == cpuc + tmp * va
        end

        uv::Tuple{T,T} = if T === BigInt
            Base.divgcd(ub, vb)
        else
            ud, vd = Base.divgcd(ub, vb)
            m = typemin(T)
            M = typemax(T)
            if !(m < ud < M && m < vd < M)
                push!(indices, j)
                return false
            end
            (ud % T, vd % T)
        end

        @inbounds Z[j] = (-2*isodd(k)+1) * Base.checked_den(uv[1], uv[2])
        # @show Z[j]
        # @assert mod((-1)^isodd(k) * ub, h) == mod(vb * x̄[j], h)
    end
    return true
end

function dixon_p(::Val{N}, A, C::Factorization{Modulo{p,T}}, Y) where {N,p,T}
    λs = [norm(x) for x in eachcol(A)]
    append!(λs, [norm(x) for x in eachcol(Y)])
    partialsort!(λs, N)
    for _ in 1:N
        popfirst!(λs)
    end
    δ::BigFloat = @static if VERSION < v"1.6-"
        (isempty(λs) ? BigFloat(1) : prod(BigFloat, λs))::BigFloat
    else
        prod(BigFloat, λs; init=one(BigFloat))
    end
    # @show δ
    # @show p
    m = ceil(Int, 2*log(δ / (MathConstants.φ - 1))/log(p))
    @assert m ≥ 1
    # @show m
    B = copy(Y)
    x̄ = BigInt.(linsolve!(C, B))
    X = copy(x̄)
    @assert A * Modulo{p,T}.(X) == B
    h = one(BigInt) # = p^i
    tmp = BigInt()
    for i in 1:m-1
        MPZ.mul_si!(h, p)
        B .= (B .- A*Integer.(X)) .÷ p
        X .= Integer.(linsolve!(C, B))
        @assert A * Modulo{p,T}.(X) == B
        # x̄ .+= h .* X
        @inbounds for j in eachindex(x̄)
            MPZ.mul!(tmp, X[j], h)
            MPZ.add!(x̄[j], tmp)
        end
    end
    MPZ.mul_si!(h, p) # h = p^m
    @assert mod.(A * x̄, h) == mod.(Y, h)
    sqh = MPZ.sqrt(h) # h = p^{m/2}
    typeofZ = Union{Matrix{Rational{Int64}},Matrix{Rational{Int128}},Matrix{Rational{BigInt}}}
    Z::typeofZ = similar(Y, Rational{Int64})
    indices = collect(reverse(eachindex(Z)))
    success = _inner_dixon_p!(indices, Z, h, x̄, sqh, tmp)
    if !success
        Z = copyuntil(first(indices), Z, Rational{Int128})
        success = _inner_dixon_p!(indices, Z, h, x̄, sqh, tmp)
        if !success
            Z = copyuntil(first(indices), Z, Rational{BigInt})
            success = _inner_dixon_p!(indices, Z, h, x̄, sqh, tmp)
            @assert success
        end
    end

    @assert eltype(Y).(A * big.(Z)) == Y
    return Z::typeofZ
    # return hcat(zeros(Rational{Int128}, N), Rational{Int128}.(x̄)')
    # Rational{Int64} is not enough for tep for instance.
end

"""
    dixon_solve(::Val{N}, A, Y) where N

Specialized solver for the linear system `A*X = Y` where `A` is a sparse symmetric
integer `n×n` matrix and `Y` is a dense integer `n×N` matrix, using Dixon's method.

Return `X` as a matrix of `Rational{Int128}`.
"""
function dixon_solve(::Val{N}, A, Y) where N
    # @show time_ns()
    typeofB = Union{
        LU{Modulo{2147483647,Int32},SparseMatrixCSC{Modulo{2147483647,Int32},Int}},
        LU{Modulo{2147483629,Int32},SparseMatrixCSC{Modulo{2147483629,Int32},Int}},
        LU{Modulo{2147483587,Int32},SparseMatrixCSC{Modulo{2147483587,Int32},Int}}
    }
    B::typeofB = rational_lu(A, false, Modulo{2147483647,Int32})
    typeofZ = Union{Matrix{Rational{Int64}},Matrix{Rational{Int128}},Matrix{Rational{BigInt}}}
    if issuccess(B)
        Z = dixon_p(Val(N), A, B, Y)'
    else
        B = rational_lu(A, false, Modulo{2147483629,Int32})
        if issuccess(B)
            Z = dixon_p(Val(N), A, B, Y)'
        else
            B = rational_lu(A, false, Modulo{2147483587,Int32})
            if issuccess(B)
                Z = dixon_p(Val(N), A, B, Y)'
            else # The probability of this being required is *extremely* low
                return rational_solve(Val(N), A, Y)
            end
        end
    end
    return hcat(zeros(eltype(Z), N), Z)::typeofZ
end


text:
export Kraskov, genentropy

"""
## k-th nearest neighbour(kNN) based

    Kraskov(; k::Int = 1, w::Int = 0) <: EntropyEstimator

Entropy estimator based on `k`-th nearest neighbor searches[^Kraskov2004].
`w` is the [Theiler window](@ref).

!!! info
    This estimator is only available for entropy estimation.
    Probabilities cannot be obtained directly.

[^Kraskov2004]: Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information. Physical review E, 69(6), 066138.
"""
Base.@kwdef struct Kraskov <: NearestNeighborEntropyEstimator
    k::Int = 1
    w::Int = 0
end

function genentropy(x::AbstractDataset{D, T}, est::Kraskov; base::Real = MathConstants.e) where {D, T}
    N = length(x)
    ρs = maximum_neighbor_distances(x, est)
    h = -digamma(est.k) + digamma(N) + log(base, ball_volume(D)) + D/N*sum(log.(base, ρs))
    return h
end


text:
module XPP
using ODE

if VERSION < v"0.4.0"
    using Docile
end

#OS-compatibility
include("os.jl")

#Definition of Core-Types
include("typedefs.jl")

#Parse ODE-Files
include("parse.jl")

#Interface for AUTO-based bifurctaion analysis
include("auto_interface.jl")

#Write model to Auto-compatible f90-file
include("auto_parser.jl")

#Routines for running simulations using XPP
include("run.jl")

#Model I/O - JSON
include("io.jl")

#Plotting Functions
include("plot.jl")

#Utility Functions
include("utils.jl")

end

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl