text:
# Julia wrapper for header: ./polylib-5.22.5/include/polylib/polylib64.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function dump_exception_stack_to_file(arg1::Ptr{FILE})
    ccall((:dump_exception_stack_to_file,polylib),Void,(Ptr{FILE},),arg1)
end

function dump_exception_stack()
    ccall((:dump_exception_stack,polylib),Void,())
end

function push_exception_on_stack(arg1::Cint,arg2::Ptr{UInt8},arg3::Ptr{UInt8},arg4::Cint)
    ccall((:push_exception_on_stack,polylib),Ptr{jmp_buf},(Cint,Ptr{UInt8},Ptr{UInt8},Cint),arg1,arg2,arg3,arg4)
end

function pop_exception_from_stack(arg1::Cint,arg2::Ptr{UInt8},arg3::Ptr{UInt8},arg4::Cint)
    ccall((:pop_exception_from_stack,polylib),Void,(Cint,Ptr{UInt8},Ptr{UInt8},Cint),arg1,arg2,arg3,arg4)
end

function throw_exception(arg1::Cint,arg2::Ptr{UInt8},arg3::Ptr{UInt8},arg4::Cint)
    ccall((:throw_exception,polylib),Void,(Cint,Ptr{UInt8},Ptr{UInt8},Cint),arg1,arg2,arg3,arg4)
end

function errormsg1(f::Ptr{UInt8},msgname::Ptr{UInt8},msg::Ptr{UInt8})
    ccall((:errormsg1,polylib),Void,(Ptr{UInt8},Ptr{UInt8},Ptr{UInt8}),f,msgname,msg)
end

function ConstraintSimplify(old::Ptr{Value},n::Ptr{Value},len::Cint,v::Ptr{Value})
    ccall((:ConstraintSimplify,polylib),Cint,(Ptr{Value},Ptr{Value},Cint,Ptr{Value}),old,n,len,v)
end

function value_alloc(want::Cint,got::Ptr{Cint})
    ccall((:value_alloc,polylib),Ptr{Value},(Cint,Ptr{Cint}),want,got)
end

function value_free(p::Ptr{Value},size::Cint)
    ccall((:value_free,polylib),Void,(Ptr{Value},Cint),p,size)
end

function Factorial(n::Cint,result::Ptr{Value})
    ccall((:Factorial,polylib),Void,(Cint,Ptr{Value}),n,result)
end

function Binomial(n::Cint,p::Cint,result::Ptr{Value})
    ccall((:Binomial,polylib),Void,(Cint,Cint,Ptr{Value}),n,p,result)
end

function CNP(a::Cint,b::Cint,result::Ptr{Value})
    ccall((:CNP,polylib),Void,(Cint,Cint,Ptr{Value}),a,b,result)
end

function Gcd(a::Value,b::Value,result::Ptr{Value})
    ccall((:Gcd,polylib),Void,(Value,Value,Ptr{Value}),a,b,result)
end

function First_Non_Zero(p::Ptr{Value},length::UInt32)
    ccall((:First_Non_Zero,polylib),Cint,(Ptr{Value},UInt32),p,length)
end

function Vector_Alloc(length::UInt32)
    ccall((:Vector_Alloc,polylib),Ptr{Vector},(UInt32,),length)
end

function Vector_Free(vector::Ptr{Vector})
    ccall((:Vector_Free,polylib),Void,(Ptr{Vector},),vector)
end

function Vector_Print(Dst::Ptr{FILE},Format::Ptr{UInt8},Vec::Ptr{Vector})
    ccall((:Vector_Print,polylib),Void,(Ptr{FILE},Ptr{UInt8},Ptr{Vector}),Dst,Format,Vec)
end

function Vector_Read()
    ccall((:Vector_Read,polylib),Ptr{Vector},())
end

function Vector_Set(p::Ptr{Value},n::Cint,length::UInt32)
    ccall((:Vector_Set,polylib),Void,(Ptr{Value},Cint,UInt32),p,n,length)
end

function Vector_Exchange(p1::Ptr{Value},p2::Ptr{Value},length::UInt32)
    ccall((:Vector_Exchange,polylib),Void,(Ptr{Value},Ptr{Value},UInt32),p1,p2,length)
end

function Vector_Copy(p1::Ptr{Value},p2::Ptr{Value},length::UInt32)
    ccall((:Vector_Copy,polylib),Void,(Ptr{Value},Ptr{Value},UInt32),p1,p2,length)
end

function Vector_Add(p1::Ptr{Value},p2::Ptr{Value},p3::Ptr{Value},length::UInt32)
    ccall((:Vector_Add,polylib),Void,(Ptr{Value},Ptr{Value},Ptr{Value},UInt32),p1,p2,p3,length)
end

function Vector_Sub(p1::Ptr{Value},p2::Ptr{Value},p3::Ptr{Value},length::UInt32)
    ccall((:Vector_Sub,polylib),Void,(Ptr{Value},Ptr{Value},Ptr{Value},UInt32),p1,p2,p3,length)
end

function Vector_Or(p1::Ptr{Value},p2::Ptr{Value},p3::Ptr{Value},length::UInt32)
    ccall((:Vector_Or,polylib),Void,(Ptr{Value},Ptr{Value},Ptr{Value},UInt32),p1,p2,p3,length)
end

function Vector_Scale(p1::Ptr{Value},p2::Ptr{Value},lambda::Value,length::UInt32)
    ccall((:Vector_Scale,polylib),Void,(Ptr{Value},Ptr{Value},Value,UInt32),p1,p2,lambda,length)
end

function Vector_AntiScale(p1::Ptr{Value},p2::Ptr{Value},lambda::Value,length::UInt32)
    ccall((:Vector_AntiScale,polylib),Void,(Ptr{Value},Ptr{Value},Value,UInt32),p1,p2,lambda,length)
end

function Vector_Oppose(p1::Ptr{Value},p2::Ptr{Value},length::UInt32)
    ccall((:Vector_Oppose,polylib),Void,(Ptr{Value},Ptr{Value},UInt32),p1,p2,length)
end

function Inner_Product(p1::Ptr{Value},p2::Ptr{Value},length::UInt32,result::Ptr{Value})
    ccall((:Inner_Product,polylib),Void,(Ptr{Value},Ptr{Value},UInt32,Ptr{Value}),p1,p2,length,result)
end

function Vector_Max(p::Ptr{Value},length::UInt32,result::Ptr{Value})
    ccall((:Vector_Max,polylib),Void,(Ptr{Value},UInt32,Ptr{Value}),p,length,result)
end

function Vector_Min(p::Ptr{Value},length::UInt32,result::Ptr{Value})
    ccall((:Vector_Min,polylib),Void,(Ptr{Value},UInt32,Ptr{Value}),p,length,result)
end

function Vector_Combine(p1::Ptr{Value},p2::Ptr{Value},p3::Ptr{Value},lambda::Value,mu::Value,length::UInt32)
    ccall((:Vector_Combine,polylib),Void,(Ptr{Value},Ptr{Value},Ptr{Value},Value,Value,UInt32),p1,p2,p3,lambda,mu,length)
end

function Vector_Equal(Vec1::Ptr{Value},Vec2::Ptr{Value},n::UInt32)
    ccall((:Vector_Equal,polylib),Cint,(Ptr{Value},Ptr{Value},UInt32),Vec1,Vec2,n)
end

function Vector_Min_Not_Zero(p::Ptr{Value},length::UInt32,index::Ptr{Cint},result::Ptr{Value})
    ccall((:Vector_Min_Not_Zero,polylib),Void,(Ptr{Value},UInt32,Ptr{Cint},Ptr{Value}),p,length,index,result)
end

function Vector_Gcd(p::Ptr{Value},length::UInt32,result::Ptr{Value})
    ccall((:Vector_Gcd,polylib),Void,(Ptr{Value},UInt32,Ptr{Value}),p,length,result)
end

function Vector_Map(p1::Ptr{Value},p2::Ptr{Value},p3::Ptr{Value},length::UInt32,f::Ptr{Void})
    ccall((:Vector_Map,polylib),Void,(Ptr{Value},Ptr{Value},Ptr{Value},UInt32,Ptr{Void}),p1,p2,p3,length,f)
end

function Vector_Normalize(p::Ptr{Value},length::UInt32)
    ccall((:Vector_Normalize,polylib),Void,(Ptr{Value},UInt32),p,length)
end

function Vector_Normalize_Positive(p::Ptr{Value},length::Cint,pos::Cint)
    ccall((:Vector_Normalize_Positive,polylib),Void,(Ptr{Value},Cint,Cint),p,length,pos)
end

function Vector_Reduce(p::Ptr{Value},length::UInt32,f::Ptr{Void},result::Ptr{Value})
    ccall((:Vector_Reduce,polylib),Void,(Ptr{Value},UInt32,Ptr{Void},Ptr{Value}),p,length,f,result)
end

function Vector_Sort(vector::Ptr{Value},n::UInt32)
    ccall((:Vector_Sort,polylib),Void,(Ptr{Value},UInt32),vector,n)
end

function Vector_IsZero(v::Ptr{Value},length::UInt32)
    ccall((:Vector_IsZero,polylib),Cint,(Ptr{Value},UInt32),v,length)
end

function Matrix_Alloc(NbRows::UInt32,NbColumns::UInt32)
    ccall((:Matrix_Alloc,polylib),Ptr{Matrix},(UInt32,UInt32),NbRows,NbColumns)
end

function Matrix_Free(Mat::Ptr{Matrix})
    ccall((:Matrix_Free,polylib),Void,(Ptr{Matrix},),Mat)
end

function Matrix_Extend(Mat::Ptr{Matrix},NbRows::UInt32)
    ccall((:Matrix_Extend,polylib),Void,(Ptr{Matrix},UInt32),Mat,NbRows)
end

function Matrix_Print(Dst::Ptr{FILE},Format::Ptr{UInt8},Mat::Ptr{Matrix})
    ccall((:Matrix_Print,polylib),Void,(Ptr{FILE},Ptr{UInt8},Ptr{Matrix}),Dst,Format,Mat)
end

function Matrix_Read_Input(Mat::Ptr{Matrix})
    ccall((:Matrix_Read_Input,polylib),Void,(Ptr{Matrix},),Mat)
end

function Matrix_Read()
    ccall((:Matrix_Read,polylib),Ptr{Matrix},())
end

function right_hermite(A::Ptr{Matrix},Hp::Ptr{Ptr{Matrix}},Up::Ptr{Ptr{Matrix}},Qp::Ptr{Ptr{Matrix}})
    ccall((:right_hermite,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),A,Hp,Up,Qp)
end

function left_hermite(A::Ptr{Matrix},Hp::Ptr{Ptr{Matrix}},Qp::Ptr{Ptr{Matrix}},Up::Ptr{Ptr{Matrix}})
    ccall((:left_hermite,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),A,Hp,Qp,Up)
end

function MatInverse(M::Ptr{Matrix},MInv::Ptr{Matrix})
    ccall((:MatInverse,polylib),Cint,(Ptr{Matrix},Ptr{Matrix}),M,MInv)
end

function rat_prodmat(S::Ptr{Matrix},X::Ptr{Matrix},P::Ptr{Matrix})
    ccall((:rat_prodmat,polylib),Void,(Ptr{Matrix},Ptr{Matrix},Ptr{Matrix}),S,X,P)
end

function Matrix_Vector_Product(mat::Ptr{Matrix},p1::Ptr{Value},p2::Ptr{Value})
    ccall((:Matrix_Vector_Product,polylib),Void,(Ptr{Matrix},Ptr{Value},Ptr{Value}),mat,p1,p2)
end

function Vector_Matrix_Product(p1::Ptr{Value},mat::Ptr{Matrix},p2::Ptr{Value})
    ccall((:Vector_Matrix_Product,polylib),Void,(Ptr{Value},Ptr{Matrix},Ptr{Value}),p1,mat,p2)
end

function Matrix_Product(mat1::Ptr{Matrix},mat2::Ptr{Matrix},mat3::Ptr{Matrix})
    ccall((:Matrix_Product,polylib),Void,(Ptr{Matrix},Ptr{Matrix},Ptr{Matrix}),mat1,mat2,mat3)
end

function Matrix_Inverse(Mat::Ptr{Matrix},MatInv::Ptr{Matrix})
    ccall((:Matrix_Inverse,polylib),Cint,(Ptr{Matrix},Ptr{Matrix}),Mat,MatInv)
end

function Polyhedron_Compute_Dual(P::Ptr{Polyhedron})
    ccall((:Polyhedron_Compute_Dual,polylib),Void,(Ptr{Polyhedron},),P)
end

function DomainConstraintSimplify(P::Ptr{Polyhedron},MaxRays::UInt32)
    ccall((:DomainConstraintSimplify,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},UInt32),P,MaxRays)
end

function AddConstraints(Con::Ptr{Value},NbConstraints::UInt32,Pol::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:AddConstraints,polylib),Ptr{Polyhedron},(Ptr{Value},UInt32,Ptr{Polyhedron},UInt32),Con,NbConstraints,Pol,NbMaxRays)
end

function AddPolyToDomain(Pol::Ptr{Polyhedron},PolDomain::Ptr{Polyhedron})
    ccall((:AddPolyToDomain,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron}),Pol,PolDomain)
end

function AddRays(Ray2::Ptr{Value},NbRay2::UInt32,Pol1::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:AddRays,polylib),Ptr{Polyhedron},(Ptr{Value},UInt32,Ptr{Polyhedron},UInt32),Ray2,NbRay2,Pol1,NbMaxRays)
end

function align_context(Pol::Ptr{Polyhedron},align_dimension::Cint,NbMaxRays::Cint)
    ccall((:align_context,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Cint,Cint),Pol,align_dimension,NbMaxRays)
end

function Constraints2Polyhedron(Constraints::Ptr{Matrix},NbMaxRays::UInt32)
    ccall((:Constraints2Polyhedron,polylib),Ptr{Polyhedron},(Ptr{Matrix},UInt32),Constraints,NbMaxRays)
end

function Disjoint_Domain(Pol::Ptr{Polyhedron},flag::Cint,NbMaxRays::UInt32)
    ccall((:Disjoint_Domain,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Cint,UInt32),Pol,flag,NbMaxRays)
end

function DomainAddConstraints(Pol::Ptr{Polyhedron},Mat::Ptr{Matrix},NbMaxRays::UInt32)
    ccall((:DomainAddConstraints,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Matrix},UInt32),Pol,Mat,NbMaxRays)
end

function DomainAddRays(Pol::Ptr{Polyhedron},Ray::Ptr{Matrix},NbMaxConstrs::UInt32)
    ccall((:DomainAddRays,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Matrix},UInt32),Pol,Ray,NbMaxConstrs)
end

function DomainConvex(Pol::Ptr{Polyhedron},NbMaxConstrs::UInt32)
    ccall((:DomainConvex,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},UInt32),Pol,NbMaxConstrs)
end

function DomainCost(Pol::Ptr{Polyhedron},Cost::Ptr{Value})
    ccall((:DomainCost,polylib),Ptr{Interval},(Ptr{Polyhedron},Ptr{Value}),Pol,Cost)
end

function DomainDifference(Pol1::Ptr{Polyhedron},Pol2::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:DomainDifference,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32),Pol1,Pol2,NbMaxRays)
end

function DomainImage(Pol::Ptr{Polyhedron},Func::Ptr{Matrix},NbMaxConstrs::UInt32)
    ccall((:DomainImage,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Matrix},UInt32),Pol,Func,NbMaxConstrs)
end

function DomainIntersection(Pol1::Ptr{Polyhedron},Pol2::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:DomainIntersection,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32),Pol1,Pol2,NbMaxRays)
end

function DomainPreimage(Pol::Ptr{Polyhedron},Func::Ptr{Matrix},NbMaxRays::UInt32)
    ccall((:DomainPreimage,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Matrix},UInt32),Pol,Func,NbMaxRays)
end

function DomainSimplify(Pol1::Ptr{Polyhedron},Pol2::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:DomainSimplify,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32),Pol1,Pol2,NbMaxRays)
end

function Stras_DomainSimplify(Pol1::Ptr{Polyhedron},Pol2::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:Stras_DomainSimplify,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32),Pol1,Pol2,NbMaxRays)
end

function DomainUnion(Pol1::Ptr{Polyhedron},Pol2::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:DomainUnion,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32),Pol1,Pol2,NbMaxRays)
end

function Domain_Copy(Pol::Ptr{Polyhedron})
    ccall((:Domain_Copy,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},),Pol)
end

function Domain_Free(Pol::Ptr{Polyhedron})
    ccall((:Domain_Free,polylib),Void,(Ptr{Polyhedron},),Pol)
end

function Domain_PrintConstraints(Dst::Ptr{FILE},Format::Ptr{UInt8},Pol::Ptr{Polyhedron})
    ccall((:Domain_PrintConstraints,polylib),Void,(Ptr{FILE},Ptr{UInt8},Ptr{Polyhedron}),Dst,Format,Pol)
end

function Empty_Polyhedron(Dimension::UInt32)
    ccall((:Empty_Polyhedron,polylib),Ptr{Polyhedron},(UInt32,),Dimension)
end

function Gauss(Mat::Ptr{Matrix},NbEq::Cint,Dimension::Cint)
    ccall((:Gauss,polylib),Cint,(Ptr{Matrix},Cint,Cint),Mat,NbEq,Dimension)
end

function lower_upper_bounds(pos::Cint,P::Ptr{Polyhedron},context::Ptr{Value},LBp::Ptr{Value},UBp::Ptr{Value})
    ccall((:lower_upper_bounds,polylib),Cint,(Cint,Ptr{Polyhedron},Ptr{Value},Ptr{Value},Ptr{Value}),pos,P,context,LBp,UBp)
end

function Polyhedron2Constraints(Pol::Ptr{Polyhedron})
    ccall((:Polyhedron2Constraints,polylib),Ptr{Matrix},(Ptr{Polyhedron},),Pol)
end

function Polyhedron2Rays(Pol::Ptr{Polyhedron})
    ccall((:Polyhedron2Rays,polylib),Ptr{Matrix},(Ptr{Polyhedron},),Pol)
end

function PolyhedronIncludes(Pol1::Ptr{Polyhedron},Pol2::Ptr{Polyhedron})
    ccall((:PolyhedronIncludes,polylib),Cint,(Ptr{Polyhedron},Ptr{Polyhedron}),Pol1,Pol2)
end

function Polyhedron_Alloc(Dimension::UInt32,NbConstraints::UInt32,NbRays::UInt32)
    ccall((:Polyhedron_Alloc,polylib),Ptr{Polyhedron},(UInt32,UInt32,UInt32),Dimension,NbConstraints,NbRays)
end

function Polyhedron_Copy(Pol::Ptr{Polyhedron})
    ccall((:Polyhedron_Copy,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},),Pol)
end

function Polyhedron_Free(Pol::Ptr{Polyhedron})
    ccall((:Polyhedron_Free,polylib),Void,(Ptr{Polyhedron},),Pol)
end

function Polyhedron_Image(Pol::Ptr{Polyhedron},Func::Ptr{Matrix},NbMaxConstrs::UInt32)
    ccall((:Polyhedron_Image,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Matrix},UInt32),Pol,Func,NbMaxConstrs)
end

function Polyhedron_Preimage(Pol::Ptr{Polyhedron},Func::Ptr{Matrix},NbMaxRays::UInt32)
    ccall((:Polyhedron_Preimage,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Matrix},UInt32),Pol,Func,NbMaxRays)
end

function Polyhedron_Print(Dst::Ptr{FILE},Format::Ptr{UInt8},Pol::Ptr{Polyhedron})
    ccall((:Polyhedron_Print,polylib),Void,(Ptr{FILE},Ptr{UInt8},Ptr{Polyhedron}),Dst,Format,Pol)
end

function Polyhedron_PrintConstraints(Dst::Ptr{FILE},Format::Ptr{UInt8},Pol::Ptr{Polyhedron})
    ccall((:Polyhedron_PrintConstraints,polylib),Void,(Ptr{FILE},Ptr{UInt8},Ptr{Polyhedron}),Dst,Format,Pol)
end

function Polyhedron_Scan(D::Ptr{Polyhedron},C::Ptr{Polyhedron},MAXRAYS::UInt32)
    ccall((:Polyhedron_Scan,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32),D,C,MAXRAYS)
end

function PolyPrint(Pol::Ptr{Polyhedron})
    ccall((:PolyPrint,polylib),Void,(Ptr{Polyhedron},),Pol)
end

function Rays2Polyhedron(Ray::Ptr{Matrix},NbMaxConstrs::UInt32)
    ccall((:Rays2Polyhedron,polylib),Ptr{Polyhedron},(Ptr{Matrix},UInt32),Ray,NbMaxConstrs)
end

function SubConstraint(Con::Ptr{Value},Pol::Ptr{Polyhedron},NbMaxRays::UInt32,Pass::Cint)
    ccall((:SubConstraint,polylib),Ptr{Polyhedron},(Ptr{Value},Ptr{Polyhedron},UInt32,Cint),Con,Pol,NbMaxRays,Pass)
end

function Universe_Polyhedron(Dimension::UInt32)
    ccall((:Universe_Polyhedron,polylib),Ptr{Polyhedron},(UInt32,),Dimension)
end

function Compute_PDomains(PD::Ptr{Param_Domain},nb_domains::Cint,working_space::Cint)
    ccall((:Compute_PDomains,polylib),Void,(Ptr{Param_Domain},Cint,Cint),PD,nb_domains,working_space)
end

function GenParamPolyhedron(Pol::Ptr{Polyhedron},Rays::Ptr{Matrix})
    ccall((:GenParamPolyhedron,polylib),Ptr{Param_Polyhedron},(Ptr{Polyhedron},Ptr{Matrix}),Pol,Rays)
end

function Param_Domain_Free(PD::Ptr{Param_Domain})
    ccall((:Param_Domain_Free,polylib),Void,(Ptr{Param_Domain},),PD)
end

function Param_Polyhedron_Free(P::Ptr{Param_Polyhedron})
    ccall((:Param_Polyhedron_Free,polylib),Void,(Ptr{Param_Polyhedron},),P)
end

function Param_Vertices_Free(PV::Ptr{Param_Vertices})
    ccall((:Param_Vertices_Free,polylib),Void,(Ptr{Param_Vertices},),PV)
end

function Param_Vertices_Print(DST::Ptr{FILE},PV::Ptr{Param_Vertices},param_names::Ptr{Ptr{UInt8}})
    ccall((:Param_Vertices_Print,polylib),Void,(Ptr{FILE},Ptr{Param_Vertices},Ptr{Ptr{UInt8}}),DST,PV,param_names)
end

function PDomainDifference(Pol1::Ptr{Polyhedron},Pol2::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:PDomainDifference,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32),Pol1,Pol2,NbMaxRays)
end

function PDomainIntersection(Pol1::Ptr{Polyhedron},Pol2::Ptr{Polyhedron},NbMaxRays::UInt32)
    ccall((:PDomainIntersection,polylib),Ptr{Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32),Pol1,Pol2,NbMaxRays)
end

function Polyhedron2Param_Domain(Din::Ptr{Polyhedron},Cin::Ptr{Polyhedron},working_space::Cint)
    ccall((:Polyhedron2Param_Domain,polylib),Ptr{Param_Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},Cint),Din,Cin,working_space)
end

function Polyhedron2Param_SimplifiedDomain(Din::Ptr{Ptr{Polyhedron}},Cin::Ptr{Polyhedron},working_space::Cint,CEq::Ptr{Ptr{Polyhedron}},CT::Ptr{Ptr{Matrix}})
    ccall((:Polyhedron2Param_SimplifiedDomain,polylib),Ptr{Param_Polyhedron},(Ptr{Ptr{Polyhedron}},Ptr{Polyhedron},Cint,Ptr{Ptr{Polyhedron}},Ptr{Ptr{Matrix}}),Din,Cin,working_space,CEq,CT)
end

function Polyhedron2Param_Vertices(Din::Ptr{Polyhedron},Cin::Ptr{Polyhedron},working_space::Cint)
    ccall((:Polyhedron2Param_Vertices,polylib),Ptr{Param_Polyhedron},(Ptr{Polyhedron},Ptr{Polyhedron},Cint),Din,Cin,working_space)
end

function Print_Domain(DST::Ptr{FILE},D::Ptr{Polyhedron},param_names::Ptr{Ptr{UInt8}})
    ccall((:Print_Domain,polylib),Void,(Ptr{FILE},Ptr{Polyhedron},Ptr{Ptr{UInt8}}),DST,D,param_names)
end

function Print_Vertex(DST::Ptr{FILE},V::Ptr{Matrix},param_names::Ptr{Ptr{UInt8}})
    ccall((:Print_Vertex,polylib),Void,(Ptr{FILE},Ptr{Matrix},Ptr{Ptr{UInt8}}),DST,V,param_names)
end

function VertexCT(V::Ptr{Matrix},CT::Ptr{Matrix})
    ccall((:VertexCT,polylib),Ptr{Matrix},(Ptr{Matrix},Ptr{Matrix}),V,CT)
end

function Param_Polyhedron_Scale_Integer(PP::Ptr{Param_Polyhedron},P::Ptr{Ptr{Polyhedron}},det::Ptr{Value},MaxRays::UInt32)
    ccall((:Param_Polyhedron_Scale_Integer,polylib),Void,(Ptr{Param_Polyhedron},Ptr{Ptr{Polyhedron}},Ptr{Value},UInt32),PP,P,det,MaxRays)
end

function Read_ParamNames(_in::Ptr{FILE},m::Cint)
    ccall((:Read_ParamNames,polylib),Ptr{Ptr{UInt8}},(Ptr{FILE},Cint),_in,m)
end

function Free_ParamNames(params::Ptr{Ptr{UInt8}},m::Cint)
    ccall((:Free_ParamNames,polylib),Void,(Ptr{Ptr{UInt8}},Cint),params,m)
end

function GaussSimplify(M::Ptr{Matrix},M2::Ptr{Matrix})
    ccall((:GaussSimplify,polylib),Cint,(Ptr{Matrix},Ptr{Matrix}),M,M2)
end

function PolyhedronLTQ(P1::Ptr{Polyhedron},P2::Ptr{Polyhedron},INDEX::Cint,PDIM::Cint,MAXRAYS::Cint)
    ccall((:PolyhedronLTQ,polylib),Cint,(Ptr{Polyhedron},Ptr{Polyhedron},Cint,Cint,Cint),P1,P2,INDEX,PDIM,MAXRAYS)
end

function PolyhedronTSort(L::Ptr{Ptr{Polyhedron}},n::UInt32,index::UInt32,pdim::UInt32,time::Ptr{Cint},pvect::Ptr{Cint},MAXRAYS::UInt32)
    ccall((:PolyhedronTSort,polylib),Cint,(Ptr{Ptr{Polyhedron}},UInt32,UInt32,UInt32,Ptr{Cint},Ptr{Cint},UInt32),L,n,index,pdim,time,pvect,MAXRAYS)
end

function Polyhedron_Not_Empty(P::Ptr{Polyhedron},C::Ptr{Polyhedron},MAXRAYS::Cint)
    ccall((:Polyhedron_Not_Empty,polylib),Cint,(Ptr{Polyhedron},Ptr{Polyhedron},Cint),P,C,MAXRAYS)
end

function count_points(pos::Cint,P::Ptr{Polyhedron},context::Ptr{Value},res::Ptr{Value})
    ccall((:count_points,polylib),Void,(Cint,Ptr{Polyhedron},Ptr{Value},Ptr{Value}),pos,P,context,res)
end

function eadd(e1::Ptr{evalue},res::Ptr{evalue})
    ccall((:eadd,polylib),Void,(Ptr{evalue},Ptr{evalue}),e1,res)
end

function ecopy(e::Ptr{enode})
    ccall((:ecopy,polylib),Ptr{enode},(Ptr{enode},),e)
end

function edot(v1::Ptr{enode},v2::Ptr{enode},res::Ptr{evalue})
    ccall((:edot,polylib),Void,(Ptr{enode},Ptr{enode},Ptr{evalue}),v1,v2,res)
end

function new_enode(_type::enode_type,size::Cint,pos::Cint)
    ccall((:new_enode,polylib),Ptr{enode},(enode_type,Cint,Cint),_type,size,pos)
end

function free_evalue_refs(e::Ptr{evalue})
    ccall((:free_evalue_refs,polylib),Void,(Ptr{evalue},),e)
end

function Polyhedron_Enumerate(P::Ptr{Polyhedron},C::Ptr{Polyhedron},MAXRAYS::UInt32,pname::Ptr{Ptr{UInt8}})
    ccall((:Polyhedron_Enumerate,polylib),Ptr{Enumeration},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32,Ptr{Ptr{UInt8}}),P,C,MAXRAYS,pname)
end

function print_enode(DST::Ptr{FILE},p::Ptr{enode},pname::Ptr{Ptr{UInt8}})
    ccall((:print_enode,polylib),Void,(Ptr{FILE},Ptr{enode},Ptr{Ptr{UInt8}}),DST,p,pname)
end

function print_evalue(DST::Ptr{FILE},e::Ptr{evalue},pname::Ptr{Ptr{UInt8}})
    ccall((:print_evalue,polylib),Void,(Ptr{FILE},Ptr{evalue},Ptr{Ptr{UInt8}}),DST,e,pname)
end

function reduce_evalue(e::Ptr{evalue})
    ccall((:reduce_evalue,polylib),Void,(Ptr{evalue},),e)
end

function Enumeration_Free(en::Ptr{Enumeration})
    ccall((:Enumeration_Free,polylib),Void,(Ptr{Enumeration},),en)
end

function Ehrhart_Quick_Apx(M::Ptr{Matrix},C::Ptr{Matrix},Validity_Lattice::Ptr{Ptr{Matrix}},MAXRAYS::UInt32)
    ccall((:Ehrhart_Quick_Apx,polylib),Ptr{Enumeration},(Ptr{Matrix},Ptr{Matrix},Ptr{Ptr{Matrix}},UInt32),M,C,Validity_Lattice,MAXRAYS)
end

function Enumeration_zero(nbParms::UInt32,maxRays::UInt32)
    ccall((:Enumeration_zero,polylib),Ptr{Enumeration},(UInt32,UInt32),nbParms,maxRays)
end

function Domain_Enumerate(D::Ptr{Polyhedron},C::Ptr{Polyhedron},MAXRAYS::UInt32,pn::Ptr{Ptr{UInt8}})
    ccall((:Domain_Enumerate,polylib),Ptr{Enumeration},(Ptr{Polyhedron},Ptr{Polyhedron},UInt32,Ptr{Ptr{UInt8}}),D,C,MAXRAYS,pn)
end

function new_eadd(e1::Ptr{evalue},res::Ptr{evalue})
    ccall((:new_eadd,polylib),Void,(Ptr{evalue},Ptr{evalue}),e1,res)
end

function compute_evalue(e::Ptr{evalue},list_args::Ptr{Value})
    ccall((:compute_evalue,polylib),Cdouble,(Ptr{evalue},Ptr{Value}),e,list_args)
end

function compute_poly(en::Ptr{Enumeration},list_args::Ptr{Value})
    ccall((:compute_poly,polylib),Ptr{Value},(Ptr{Enumeration},Ptr{Value}),en,list_args)
end

function in_domain(P::Ptr{Polyhedron},list_args::Ptr{Value})
    ccall((:in_domain,polylib),Cint,(Ptr{Polyhedron},Ptr{Value}),P,list_args)
end

function SolveDiophantine(M::Ptr{Matrix},U::Ptr{Ptr{Matrix}},X::Ptr{Ptr{Vector}})
    ccall((:SolveDiophantine,polylib),Cint,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Vector}}),M,U,X)
end

function AffineHermite(A::Ptr{Lattice},H::Ptr{Ptr{Lattice}},U::Ptr{Ptr{Matrix}})
    ccall((:AffineHermite,polylib),Void,(Ptr{Lattice},Ptr{Ptr{Lattice}},Ptr{Ptr{Matrix}}),A,H,U)
end

function AffineSmith(A::Ptr{Lattice},U::Ptr{Ptr{Lattice}},V::Ptr{Ptr{Lattice}},Diag::Ptr{Ptr{Lattice}})
    ccall((:AffineSmith,polylib),Void,(Ptr{Lattice},Ptr{Ptr{Lattice}},Ptr{Ptr{Lattice}},Ptr{Ptr{Lattice}}),A,U,V,Diag)
end

function ChangeLatticeDimension(A::Ptr{Lattice},dimension::Cint)
    ccall((:ChangeLatticeDimension,polylib),Ptr{Lattice},(Ptr{Lattice},Cint),A,dimension)
end

function EmptyLattice(dimension::Cint)
    ccall((:EmptyLattice,polylib),Ptr{Lattice},(Cint,),dimension)
end

function ExtractLinearPart(A::Ptr{Lattice})
    ccall((:ExtractLinearPart,polylib),Ptr{Lattice},(Ptr{Lattice},),A)
end

function FindHermiteBasisofDomain(A::Ptr{Polyhedron},B::Ptr{Ptr{Matrix}})
    ccall((:FindHermiteBasisofDomain,polylib),Cint,(Ptr{Polyhedron},Ptr{Ptr{Matrix}}),A,B)
end

function Homogenise(A::Ptr{Lattice},Forward::Bool)
    ccall((:Homogenise,polylib),Ptr{Lattice},(Ptr{Lattice},Bool),A,Forward)
end

function intcompare(a::Ptr{Void},b::Ptr{Void})
    ccall((:intcompare,polylib),Cint,(Ptr{Void},Ptr{Void}),a,b)
end

function isEmptyLattice(A::Ptr{Lattice})
    ccall((:isEmptyLattice,polylib),Bool,(Ptr{Lattice},),A)
end

function isfulldim(m::Ptr{Matrix})
    ccall((:isfulldim,polylib),Bool,(Ptr{Matrix},),m)
end

function IsLattice(m::Ptr{Matrix})
    ccall((:IsLattice,polylib),Bool,(Ptr{Matrix},),m)
end

function isLinear(A::Ptr{Lattice})
    ccall((:isLinear,polylib),Bool,(Ptr{Lattice},),A)
end

function LatticeDifference(A::Ptr{Lattice},B::Ptr{Lattice})
    ccall((:LatticeDifference,polylib),Ptr{LatticeUnion},(Ptr{Lattice},Ptr{Lattice}),A,B)
end

function LatticeImage(A::Ptr{Lattice},M::Ptr{Matrix})
    ccall((:LatticeImage,polylib),Ptr{Lattice},(Ptr{Lattice},Ptr{Matrix}),A,M)
end

function LatticeIncludes(A::Ptr{Lattice},B::Ptr{Lattice})
    ccall((:LatticeIncludes,polylib),Bool,(Ptr{Lattice},Ptr{Lattice}),A,B)
end

function LatticeIntersection(X::Ptr{Lattice},Y::Ptr{Lattice})
    ccall((:LatticeIntersection,polylib),Ptr{Lattice},(Ptr{Lattice},Ptr{Lattice}),X,Y)
end

function LatticePreimage(L::Ptr{Lattice},G::Ptr{Matrix})
    ccall((:LatticePreimage,polylib),Ptr{Lattice},(Ptr{Lattice},Ptr{Matrix}),L,G)
end

function LatticeSimplify(latlist::Ptr{LatticeUnion})
    ccall((:LatticeSimplify,polylib),Ptr{LatticeUnion},(Ptr{LatticeUnion},),latlist)
end

function LatticeUnion_Alloc()
    ccall((:LatticeUnion_Alloc,polylib),Ptr{LatticeUnion},())
end

function LatticeUnion_Free(Head::Ptr{LatticeUnion})
    ccall((:LatticeUnion_Free,polylib),Void,(Ptr{LatticeUnion},),Head)
end

function PrintLatticeUnion(fp::Ptr{FILE},format::Ptr{UInt8},Head::Ptr{LatticeUnion})
    ccall((:PrintLatticeUnion,polylib),Void,(Ptr{FILE},Ptr{UInt8},Ptr{LatticeUnion}),fp,format,Head)
end

function sameAffinepart(A::Ptr{Lattice},B::Ptr{Lattice})
    ccall((:sameAffinepart,polylib),Bool,(Ptr{Lattice},Ptr{Lattice}),A,B)
end

function sameLattice(A::Ptr{Lattice},B::Ptr{Lattice})
    ccall((:sameLattice,polylib),Bool,(Ptr{Lattice},Ptr{Lattice}),A,B)
end

function Lattice2LatticeUnion(X::Ptr{Lattice},Y::Ptr{Lattice})
    ccall((:Lattice2LatticeUnion,polylib),Ptr{LatticeUnion},(Ptr{Lattice},Ptr{Lattice}),X,Y)
end

function Lcm3(a::Value,b::Value,c::Ptr{Value})
    ccall((:Lcm3,polylib),Void,(Value,Value,Ptr{Value}),a,b,c)
end

function AddANullColumn(M::Ptr{Matrix})
    ccall((:AddANullColumn,polylib),Ptr{Matrix},(Ptr{Matrix},),M)
end

function AddANullRow(M::Ptr{Matrix})
    ccall((:AddANullRow,polylib),Ptr{Matrix},(Ptr{Matrix},),M)
end

function ExchangeColumns(M::Ptr{Matrix},Column1::Cint,Column2::Cint)
    ccall((:ExchangeColumns,polylib),Void,(Ptr{Matrix},Cint,Cint),M,Column1,Column2)
end

function ExchangeRows(M::Ptr{Matrix},Row1::Cint,Row2::Cint)
    ccall((:ExchangeRows,polylib),Void,(Ptr{Matrix},Cint,Cint),M,Row1,Row2)
end

function findHermiteBasis(M::Ptr{Matrix},Result::Ptr{Ptr{Matrix}})
    ccall((:findHermiteBasis,polylib),Cint,(Ptr{Matrix},Ptr{Ptr{Matrix}}),M,Result)
end

function Identity(size::UInt32)
    ccall((:Identity,polylib),Ptr{Matrix},(UInt32,),size)
end

function isinHnf(A::Ptr{Matrix})
    ccall((:isinHnf,polylib),Bool,(Ptr{Matrix},),A)
end

function isIntegral(A::Ptr{Matrix})
    ccall((:isIntegral,polylib),Bool,(Ptr{Matrix},),A)
end

function Lcm(i::Value,j::Value)
    ccall((:Lcm,polylib),Ptr{Value},(Value,Value),i,j)
end

function Matrix_Copy(Src::Ptr{Matrix})
    ccall((:Matrix_Copy,polylib),Ptr{Matrix},(Ptr{Matrix},),Src)
end

function PutColumnFirst(X::Ptr{Matrix},Columnnumber::Cint)
    ccall((:PutColumnFirst,polylib),Void,(Ptr{Matrix},Cint),X,Columnnumber)
end

function PutColumnLast(X::Ptr{Matrix},Columnnumber::Cint)
    ccall((:PutColumnLast,polylib),Void,(Ptr{Matrix},Cint),X,Columnnumber)
end

function PutRowFirst(X::Ptr{Matrix},Rownumber::Cint)
    ccall((:PutRowFirst,polylib),Void,(Ptr{Matrix},Cint),X,Rownumber)
end

function PutRowLast(X::Ptr{Matrix},Rownumber::Cint)
    ccall((:PutRowLast,polylib),Void,(Ptr{Matrix},Cint),X,Rownumber)
end

function RemoveNColumns(M::Ptr{Matrix},FirstColumnnumber::Cint,NumColumns::Cint)
    ccall((:RemoveNColumns,polylib),Ptr{Matrix},(Ptr{Matrix},Cint,Cint),M,FirstColumnnumber,NumColumns)
end

function RemoveColumn(M::Ptr{Matrix},Columnnumber::Cint)
    ccall((:RemoveColumn,polylib),Ptr{Matrix},(Ptr{Matrix},Cint),M,Columnnumber)
end

function RemoveRow(M::Ptr{Matrix},Rownumber::Cint)
    ccall((:RemoveRow,polylib),Ptr{Matrix},(Ptr{Matrix},Cint),M,Rownumber)
end

function Transpose(A::Ptr{Matrix})
    ccall((:Transpose,polylib),Ptr{Matrix},(Ptr{Matrix},),A)
end

function Smith(A::Ptr{Matrix},U::Ptr{Ptr{Matrix}},V::Ptr{Ptr{Matrix}},Product::Ptr{Ptr{Matrix}})
    ccall((:Smith,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),A,U,V,Product)
end

function Hermite(A::Ptr{Matrix},H::Ptr{Ptr{Matrix}},U::Ptr{Ptr{Matrix}})
    ccall((:Hermite,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),A,H,U)
end

function CanonicalForm(Zpol::Ptr{ZPolyhedron},Result::Ptr{Ptr{ZPolyhedron}},Basis::Ptr{Ptr{Matrix}})
    ccall((:CanonicalForm,polylib),Void,(Ptr{ZPolyhedron},Ptr{Ptr{ZPolyhedron}},Ptr{Ptr{Matrix}}),Zpol,Result,Basis)
end

function EmptyZPolyhedron(dimension::Cint)
    ccall((:EmptyZPolyhedron,polylib),Ptr{ZPolyhedron},(Cint,),dimension)
end

function IntegraliseLattice(A::Ptr{ZPolyhedron})
    ccall((:IntegraliseLattice,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},),A)
end

function isEmptyZPolyhedron(Zpol::Ptr{ZPolyhedron})
    ccall((:isEmptyZPolyhedron,polylib),Bool,(Ptr{ZPolyhedron},),Zpol)
end

function ZDomainDifference(A::Ptr{ZPolyhedron},B::Ptr{ZPolyhedron})
    ccall((:ZDomainDifference,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},Ptr{ZPolyhedron}),A,B)
end

function ZDomainImage(A::Ptr{ZPolyhedron},Func::Ptr{Matrix})
    ccall((:ZDomainImage,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},Ptr{Matrix}),A,Func)
end

function ZDomainIncludes(A::Ptr{ZPolyhedron},B::Ptr{ZPolyhedron})
    ccall((:ZDomainIncludes,polylib),Bool,(Ptr{ZPolyhedron},Ptr{ZPolyhedron}),A,B)
end

function ZDomainIntersection(A::Ptr{ZPolyhedron},B::Ptr{ZPolyhedron})
    ccall((:ZDomainIntersection,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},Ptr{ZPolyhedron}),A,B)
end

function ZDomainPreimage(A::Ptr{ZPolyhedron},Func::Ptr{Matrix})
    ccall((:ZDomainPreimage,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},Ptr{Matrix}),A,Func)
end

function ZDomainPrint(fp::Ptr{FILE},format::Ptr{UInt8},A::Ptr{ZPolyhedron})
    ccall((:ZDomainPrint,polylib),Void,(Ptr{FILE},Ptr{UInt8},Ptr{ZPolyhedron}),fp,format,A)
end

function ZDomainSimplify(ZDom::Ptr{ZPolyhedron})
    ccall((:ZDomainSimplify,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},),ZDom)
end

function ZDomainUnion(A::Ptr{ZPolyhedron},B::Ptr{ZPolyhedron})
    ccall((:ZDomainUnion,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},Ptr{ZPolyhedron}),A,B)
end

function ZDomain_Copy(Head::Ptr{ZPolyhedron})
    ccall((:ZDomain_Copy,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},),Head)
end

function ZDomain_Free(Head::Ptr{ZPolyhedron})
    ccall((:ZDomain_Free,polylib),Void,(Ptr{ZPolyhedron},),Head)
end

function ZPolyhedronIncludes(A::Ptr{ZPolyhedron},B::Ptr{ZPolyhedron})
    ccall((:ZPolyhedronIncludes,polylib),Bool,(Ptr{ZPolyhedron},Ptr{ZPolyhedron}),A,B)
end

function ZPolyhedron_Alloc(Lat::Ptr{Lattice},Poly::Ptr{Polyhedron})
    ccall((:ZPolyhedron_Alloc,polylib),Ptr{ZPolyhedron},(Ptr{Lattice},Ptr{Polyhedron}),Lat,Poly)
end

function SplitZpolyhedron(ZPol::Ptr{ZPolyhedron},B::Ptr{Lattice})
    ccall((:SplitZpolyhedron,polylib),Ptr{ZPolyhedron},(Ptr{ZPolyhedron},Ptr{Lattice}),ZPol,B)
end

function Smith(A::Ptr{Matrix},U::Ptr{Ptr{Matrix}},V::Ptr{Ptr{Matrix}},Product::Ptr{Ptr{Matrix}})
    ccall((:Smith,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),A,U,V,Product)
end

function Hermite(A::Ptr{Matrix},H::Ptr{Ptr{Matrix}},U::Ptr{Ptr{Matrix}})
    ccall((:Hermite,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),A,H,U)
end

function Smith(A::Ptr{Matrix},U::Ptr{Ptr{Matrix}},V::Ptr{Ptr{Matrix}},Product::Ptr{Ptr{Matrix}})
    ccall((:Smith,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),A,U,V,Product)
end

function Hermite(A::Ptr{Matrix},H::Ptr{Ptr{Matrix}},U::Ptr{Ptr{Matrix}})
    ccall((:Hermite,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),A,H,U)
end

function Equalities_integerSolution(Eqs::Ptr{Matrix},sol::Ptr{Ptr{Matrix}})
    ccall((:Equalities_integerSolution,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}}),Eqs,sol)
end

function Equalities_validityLattice(Eqs::Ptr{Matrix},a::Cint,vl::Ptr{Ptr{Matrix}})
    ccall((:Equalities_validityLattice,polylib),Void,(Ptr{Matrix},Cint,Ptr{Ptr{Matrix}}),Eqs,a,vl)
end

function Equalities_intModBasis(B::Ptr{Matrix},C::Ptr{Matrix},d::Ptr{Matrix},imb::Ptr{Ptr{Matrix}})
    ccall((:Equalities_intModBasis,polylib),Void,(Ptr{Matrix},Ptr{Matrix},Ptr{Matrix},Ptr{Ptr{Matrix}}),B,C,d,imb)
end

function Constraints_fullDimensionize(M::Ptr{Ptr{Matrix}},C::Ptr{Ptr{Matrix}},VL::Ptr{Ptr{Matrix}},Eqs::Ptr{Ptr{Matrix}},ParmEqs::Ptr{Ptr{Matrix}},elimVars::Ptr{Ptr{UInt32}},elimParms::Ptr{Ptr{UInt32}},maxRays::Cint)
    ccall((:Constraints_fullDimensionize,polylib),Void,(Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Ptr{Ptr{UInt32}},Ptr{Ptr{UInt32}},Cint),M,C,VL,Eqs,ParmEqs,elimVars,elimParms,maxRays)
end

function Constraints_Remove_parm_eqs(M::Ptr{Ptr{Matrix}},Ctxt::Ptr{Ptr{Matrix}},renderSpace::Cint,elimParms::Ptr{Ptr{UInt32}})
    ccall((:Constraints_Remove_parm_eqs,polylib),Ptr{Matrix},(Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}},Cint,Ptr{Ptr{UInt32}}),M,Ctxt,renderSpace,elimParms)
end

function Constraints_removeElimCols(M::Ptr{Matrix},nbVars::UInt32,elimParms::Ptr{UInt32},newM::Ptr{Ptr{Matrix}})
    ccall((:Constraints_removeElimCols,polylib),Void,(Ptr{Matrix},UInt32,Ptr{UInt32},Ptr{Ptr{Matrix}}),M,nbVars,elimParms,newM)
end

function Lattice_extractSubLattice(lat::Ptr{Matrix},k::UInt32,subLat::Ptr{Ptr{Matrix}})
    ccall((:Lattice_extractSubLattice,polylib),Void,(Ptr{Matrix},UInt32,Ptr{Ptr{Matrix}}),lat,k,subLat)
end

function Polyhedron_Remove_parm_eqs(P::Ptr{Ptr{Polyhedron}},C::Ptr{Ptr{Polyhedron}},renderSpace::Cint,elimParms::Ptr{Ptr{UInt32}},maxRays::Cint)
    ccall((:Polyhedron_Remove_parm_eqs,polylib),Ptr{Polyhedron},(Ptr{Ptr{Polyhedron}},Ptr{Ptr{Polyhedron}},Cint,Ptr{Ptr{UInt32}},Cint),P,C,renderSpace,elimParms,maxRays)
end

function int_ker(M::Ptr{Matrix})
    ccall((:int_ker,polylib),Ptr{Matrix},(Ptr{Matrix},),M)
end

function full_dimensionize(M::Ptr{Matrix},nb_parms::Cint,Validity_Lattice::Ptr{Ptr{Matrix}})
    ccall((:full_dimensionize,polylib),Ptr{Matrix},(Ptr{Matrix},Cint,Ptr{Ptr{Matrix}}),M,nb_parms,Validity_Lattice)
end

function affine_periods(M::Ptr{Matrix},d::Ptr{Matrix})
    ccall((:affine_periods,polylib),Ptr{Matrix},(Ptr{Matrix},Ptr{Matrix}),M,d)
end

function int_mod_basis(Bp::Ptr{Matrix},Cp::Ptr{Matrix},d::Ptr{Matrix})
    ccall((:int_mod_basis,polylib),Ptr{Matrix},(Ptr{Matrix},Ptr{Matrix},Ptr{Matrix}),Bp,Cp,d)
end

function compress_parms(E::Ptr{Matrix},nb_parms::Cint)
    ccall((:compress_parms,polylib),Ptr{Matrix},(Ptr{Matrix},Cint),E,nb_parms)
end

function nb_bits(x::Culonglong)
    ccall((:nb_bits,polylib),UInt32,(Culonglong,),x)
end

function permutation_inverse(perm::Ptr{UInt32},nb_elems::UInt32)
    ccall((:permutation_inverse,polylib),Ptr{UInt32},(Ptr{UInt32},UInt32),perm,nb_elems)
end

function mtransformation_permute(transf::Ptr{Matrix},permutation::Ptr{UInt32})
    ccall((:mtransformation_permute,polylib),Ptr{Matrix},(Ptr{Matrix},Ptr{UInt32}),transf,permutation)
end

function mpolyhedron_permute(polyh::Ptr{Matrix},permutation::Ptr{UInt32})
    ccall((:mpolyhedron_permute,polylib),Ptr{Matrix},(Ptr{Matrix},Ptr{UInt32}),polyh,permutation)
end

function Constraints_permute(C::Ptr{Matrix},perm::Ptr{UInt32},Cp::Ptr{Ptr{Matrix}})
    ccall((:Constraints_permute,polylib),Void,(Ptr{Matrix},Ptr{UInt32},Ptr{Ptr{Matrix}}),C,perm,Cp)
end

function eliminable_vars(Eqs::Ptr{Matrix},start::UInt32,_end::UInt32)
    ccall((:eliminable_vars,polylib),Culonglong,(Ptr{Matrix},UInt32,UInt32),Eqs,start,_end)
end

function find_a_permutation(Eqs::Ptr{Matrix},nb_parms::UInt32)
    ccall((:find_a_permutation,polylib),Ptr{UInt32},(Ptr{Matrix},UInt32),Eqs,nb_parms)
end

function permutation_for_full_dim2(vars_to_keep::Ptr{UInt32},nb_keep::UInt32,nb_vars_parms::UInt32,nb_parms::UInt32)
    ccall((:permutation_for_full_dim2,polylib),Ptr{UInt32},(Ptr{UInt32},UInt32,UInt32,UInt32),vars_to_keep,nb_keep,nb_vars_parms,nb_parms)
end

function constraintsView(P::Ptr{Polyhedron})
    ccall((:constraintsView,polylib),Ptr{Matrix},(Ptr{Polyhedron},),P)
end

function constraintsView_Free(M::Ptr{Matrix})
    ccall((:constraintsView_Free,polylib),Void,(Ptr{Matrix},),M)
end

function split_constraints(M::Ptr{Matrix},Eqs::Ptr{Ptr{Matrix}},Ineqs::Ptr{Ptr{Matrix}})
    ccall((:split_constraints,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Ptr{Matrix}}),M,Eqs,Ineqs)
end

function Identity_Matrix(dim::UInt32)
    ccall((:Identity_Matrix,polylib),Ptr{Matrix},(UInt32,),dim)
end

function Matrix_identity(dim::UInt32,I::Ptr{Ptr{Matrix}})
    ccall((:Matrix_identity,polylib),Void,(UInt32,Ptr{Ptr{Matrix}}),dim,I)
end

function mtransformation_inverse(transf::Ptr{Matrix},inv::Ptr{Ptr{Matrix}},g::Ptr{Value})
    ccall((:mtransformation_inverse,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}},Ptr{Value}),transf,inv,g)
end

function mpolyhedron_simplify(polyh::Ptr{Matrix})
    ccall((:mpolyhedron_simplify,polylib),Void,(Ptr{Matrix},),polyh)
end

function mpolyhedron_inflate(polyh::Ptr{Matrix},nb_parms::UInt32)
    ccall((:mpolyhedron_inflate,polylib),Void,(Ptr{Matrix},UInt32),polyh,nb_parms)
end

function mpolyhedron_deflate(polyh::Ptr{Matrix},nb_parms::UInt32)
    ccall((:mpolyhedron_deflate,polylib),Void,(Ptr{Matrix},UInt32),polyh,nb_parms)
end

function eliminate_var_with_constr(Eliminator::Ptr{Matrix},eliminator_row::UInt32,Victim::Ptr{Matrix},victim_row::UInt32,var_to_elim::UInt32)
    ccall((:eliminate_var_with_constr,polylib),Void,(Ptr{Matrix},UInt32,Ptr{Matrix},UInt32,UInt32),Eliminator,eliminator_row,Victim,victim_row,var_to_elim)
end

function mpolyhedron_compress_last_vars(M::Ptr{Matrix},compression::Ptr{Matrix})
    ccall((:mpolyhedron_compress_last_vars,polylib),Void,(Ptr{Matrix},Ptr{Matrix}),M,compression)
end

function mpolyhedron_eliminate_first_variables(Eqs::Ptr{Matrix},Ineqs::Ptr{Matrix})
    ccall((:mpolyhedron_eliminate_first_variables,polylib),UInt32,(Ptr{Matrix},Ptr{Matrix}),Eqs,Ineqs)
end

function Matrix_subMatrix(M::Ptr{Matrix},sr::UInt32,sc::UInt32,nbR::UInt32,nbC::UInt32,sub::Ptr{Ptr{Matrix}})
    ccall((:Matrix_subMatrix,polylib),Void,(Ptr{Matrix},UInt32,UInt32,UInt32,UInt32,Ptr{Ptr{Matrix}}),M,sr,sc,nbR,nbC,sub)
end

function Matrix_clone(M::Ptr{Matrix},Cl::Ptr{Ptr{Matrix}})
    ccall((:Matrix_clone,polylib),Void,(Ptr{Matrix},Ptr{Ptr{Matrix}}),M,Cl)
end

function Matrix_copySubMatrix(M1::Ptr{Matrix},sr1::UInt32,sc1::UInt32,nbR::UInt32,nbC::UInt32,M2::Ptr{Matrix},sr2::UInt32,sc2::UInt32)
    ccall((:Matrix_copySubMatrix,polylib),Void,(Ptr{Matrix},UInt32,UInt32,UInt32,UInt32,Ptr{Matrix},UInt32,UInt32),M1,sr1,sc1,nbR,nbC,M2,sr2,sc2)
end

function Matrix_oppose(M::Ptr{Matrix})
    ccall((:Matrix_oppose,polylib),Void,(Ptr{Matrix},),M)
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

# load pkgs
using Unitful
using MultiStateSystems

# pkg const
const _MSS = MultiStateSystems

# setting for a specific analysis
cls = SteadyStateProcess()

# include the stochastic models for the elements
stdᵃᵈ = include(joinpath(_MSS.BASE_DIR,"examples/lvdc/elements/acdc.jl"))
stdᵃᶜ = include(joinpath(_MSS.BASE_DIR,"examples/lvdc/elements/ac.jl"))
stdᵈᵈ = include(joinpath(_MSS.BASE_DIR,"examples/lvdc/elements/dcdc.jl"))
stdᵇ  = include(joinpath(_MSS.BASE_DIR,"examples/lvdc/elements/battery.jl"))

# solve the state-transitions diagrams
solve!(stdᵃᵈ, cls)
solve!(stdᵃᶜ, cls)
solve!(stdᵈᵈ, cls)
solve!(stdᵇ, cls)

# initialize the network corresponding to the low voltage dc system
ntwˡᵛ = Network()

# add the users to the ntw
add_users!(ntwˡᵛ, node = [3,4])
add_sources!(ntwˡᵛ, node = [1,3],
                    name = ["ac grid", "battery"],
                    std  = [stdᵃᶜ,stdᵇ])
add_components!(ntwˡᵛ, edge = [(1,2),(2,3),(2,4)],
                       name = ["ac/dc", "dc/dc", "dc/ac"],
                       std  = [stdᵃᵈ, stdᵈᵈ, stdᵃᵈ])

# solve the network
solve!(ntwˡᵛ)

ext: jl


text:
"""
    galilean_trajectory_sample!(m, e, tuner, cache)

Given 'm <: GMC_NS_Model', 'e <: GMC_NS_Ensemble', 'tuner <: GMC_Tuner', and 'cache <: Union{Nothing,GMC_NS_Model}', return the next model-sample from a Gailean trajectory through parameter space, as well as any cached reflection model, obeying 'new_m.log_Li >= e.contour'.

Implements Skilling's 2019 update to the GMC algorithm. Briefly: proceed "forward" along the trajectory's velocity vector unless the likelihood contour is hit; in this case, attempt to reflect "east", then "west", then "south". If a reflection is found, return the model in the original position with the new velocity vector (maintaining detailed balance), as well as the cached reflection model "forward" along this new vector. If no model with likelihood greater than contour can be found, return the original model (as long as the trajectory is alive it will be retried at lower timestep τ).

Skilling, John. “Galilean and Hamiltonian Monte Carlo.” In Proceedings, 33:8. Garching, Germany: MDPI, 2019.
"""
function galilean_trajectory_sample!(m, e, tuner, cache)
    t=m.trajectory; next_i=m.i+1; cache_i=m.i+2

    e.GMC_timestep_η > 0. ? (τ=max(tuner.τ+rand(Normal(0,tuner.τ*e.GMC_timestep_η)),eps())) : (τ=tuner.τ) #perturb τ if specified
    d=τ*m.v #distance vector for given timestep
    fwd_pos,adj_d,box_rflct=box_move(m.pos,d,e.box)

    if !box_rflct #if we're not stopped by the sampling box
        fwd_m=e.model_initλ(t, next_i, to_prior.(fwd_pos,e.priors), fwd_pos, m.v, e.obs, e.constants...)  #try to proceed along distance vector
        if m.θ==fwd_m.θ #if forward motion is no longer making a difference to the parameter vector (ie timestep is too small), kill the trajectory and bail out
            tuner.τ=e.GMC_τ_death
            return m, nothing
        end
    else
        fwd_m=m #placeholder, no effect on box reflection
    end

    if box_rflct || (fwd_m.log_Li < m.log_Li) #if no model is available from the position along the distance vector, search for reflections and store in cache
        process_report!(tuner, false) #report a reflection off a lh contour

        lhδ=lps(m.log_Li,-fwd_m.log_Li) #get the likelihood difference between the two points
        n=boundary_norm(adj_d,lhδ) #get the gradient normal for the distance and lhδ
        box_rflct ? (v′=box_reflect(m.pos,e.box,m.v,e.GMC_reflect_η)) : v′=reflect(m.v,n,e.GMC_reflect_η) #get the reflected velocity vector off the boundary "east", or off the sampling box
        rd=τ*v′ #reflected distance given the timestep
        east_pos,_=box_move(m.pos,rd,e.box)
        cache=e.model_initλ(t, cache_i, to_prior.(east_pos,e.priors), east_pos, v′, e.obs, e.constants...) #try going>

        if cache.log_Li < m.log_Li && !box_rflct #no west reflection off the sampling box
            process_report!(tuner, false)

            west_pos,_=box_move(m.pos,-rd,e.box)
            cache=e.model_initλ(t, cache_i, to_prior.(west_pos,e.priors), west_pos, -v′, e.obs, e.constants...)  #if east reflection fails, try west
        end

        if cache.log_Li < m.log_Li
            process_report!(tuner, false)

            e.GMC_reflect_η > 0. ? (v′=r_perturb(-m.v, e.GMC_reflect_η)) : (v′=-m.v)
            sd=τ*v′ #south distance along perturbed reflexn vec (inverted m.v)

            south_pos,_=box_move(m.pos,sd,e.box)
            cache=e.model_initλ(t, cache_i, to_prior.(south_pos,e.priors),south_pos,v′, e.obs, e.constants...)  #if west or box reflection fails, try south (reversed along the particle's vector)
        end
        cache.log_Li < m.log_Li && (process_report!(tuner, false)) #if that fails wait for smaller τ
    end

    if cache!==nothing && cache.log_Li >= m.log_Li #a reflection that results in a model above contour exists, so return a model in place with the appropriate velocity vector and the cached model down the reflection vector
        r_m=construct_reflection(m, next_i, cache.v)
        return (r_m, cache)
    else #if no reflections occurred or no reflected positions above the contour were found - either forward search succeeded or the whole search failed
        cache!==nothing && (cache=nothing) #reset cache if it's been set with a model below contour
        fwd_m.i == next_i && fwd_m.log_Li >= m.log_Li ? (process_report!(tuner, true); return fwd_m, cache) : (return m, cache) #if the forward model is a new iterate and at or above the contour, return that, otherwise return the original model
    end
end

                function box_move(pos,d,box)
                    if !(all(box[:,1].<pos+d.<box[:,2])) #if not all coordinates after move would be in the box
                        boundary_d=box.-pos #get distances to box boundaries
                        b=isapprox.(pos,box) #find any boundaries that have been exceeded due to numerical creep
                        boundary_d[b].=0 #correct these distances
                        boundary_t=boundary_d./d #get time to hit boundaries, given particle displacment d

                        for idx in findall(x->x>=0.,boundary_t) #for any boundaries the particle is riding (distance 0), give a small positive time if the particle is headed into the boundary to allow that dimension to be selected for bounding
                            (dim,bound)=(idx[1],idx[2])
                            boundary_t[idx]==0. && ((bound==1 && d[dim] < 0.) || (bound==2 && d[dim] > 0.)) && (boundary_t[idx]=nextfloat(0.))
                        end

                        first_b_idx=findfirst(isequal(minimum(boundary_t[findall(x->x>0.,boundary_t)])),boundary_t) #find the index of the first boundary that would be hit
                        first_b_d=boundary_d[first_b_idx]
                        d=d.*(first_b_d/d[first_b_idx[1]])
                    end

                    all(d.==0.) ? box_reflect=true : box_reflect=false

                    new_pos=pos+d
                    !(all(box[:,1].<new_pos.<box[:,2])) && box_bound!(new_pos,box)

                    return new_pos, d, box_reflect
                end

                function box_reflect(pos,box,v,η)
                    b=isapprox.(pos,box)
                    low_idxs=[v[i]<0&&b[i,1] for i in 1:length(v)]
                    hi_idxs=[v[i]>0&&b[i,2] for i in 1:length(v)]

                    v′=copy(v)
                    v′[low_idxs]=-v[low_idxs]
                    v′[hi_idxs]=-v[hi_idxs]

                    η > 0. ? (return r_perturb(v′,η)) : (return v′)
                end

                function boundary_norm(v, lhδ)
                    b=lhδ./v
                    b[b.==Inf].=prevfloat(Inf) #rectify overflow
                    b[b.==0.].=rand([nextfloat(0.),prevfloat(0.)]) #rectify exact 0.
                    b[b.==-Inf].=nextfloat(-Inf) #rectify underflow
                    return normalize(b)
                end

                function reflect(v, n, η)
                    v′=v-(2*n*dot(n,v))
                    η > 0. ? (return r_perturb(v′,η)) : (return v′)
                end

                function r_perturb(v′, η)
                    return v′*cos(η) + rand(MvNormal(length(v′),1.))*sin(η)
                end


text:
module LASArrays

using FileIO, MappedArrays, StaticArrays, LinearAlgebra, CoordinateTransformations

abstract type AbstractPointDataRecordFormat end

struct PointDataRecordFormat0 <: AbstractPointDataRecordFormat
    x::Int32
    y::Int32
    z::Int32
    intensity::UInt16
    flag::UInt8
    classification::UInt8
    scanangle::Int8
    userdata::UInt8
    ptsrcid::UInt16
end

struct PointDataRecordFormat1 <: AbstractPointDataRecordFormat
    x::Int32
    y::Int32
    z::Int32
    intensity::UInt16
    flag::UInt8
    classification::UInt8
    scanangle::Int8
    userdata::UInt8
    ptsrcid::UInt16
    gpstime::Float64
end

struct PointDataRecordFormat2 <: AbstractPointDataRecordFormat
    x::Int32
    y::Int32
    z::Int32
    intensity::UInt16
    flag::UInt8
    classification::UInt8
    scanangle::Int8
    userdata::UInt8
    ptsrcid::UInt16
    red::UInt16
    green::UInt16
    blue::UInt16
end

struct PointDataRecordFormat3 <: AbstractPointDataRecordFormat
    x::Int32
    y::Int32
    z::Int32
    intensity::UInt16
    flag::UInt8
    classification::UInt8
    scanangle::Int8
    userdata::UInt8
    ptsrcid::UInt16
    gpstime::Float64
    red::UInt16
    green::UInt16
    blue::UInt16
end

#header offsets - number of bytes from beginning of file for a given header field including magic bytes(LASF)
const POINTDATA_OFFSET_SEEK = 96
const N_VLR_SEEK = 100
const FORMAT_SEEK = 104
const N_RECORDS_SEEK = 105
const X_SCALE_SEEK = 131
#probably don't need these since they are a contiguous block, but maybe useful later?
#const Y_SCALE_SEEK = 139
#const Z_SCALE_SEEK = 147
#const X_OFFSET_SEEK = 155
#const Y_OFFSET_SEEK = 163
#const Z_OFFSET_SEEK = 171

#point data record offsets - number of bytes from beginning of file to start of point records
pointdata_start(io::Stream{format"LAS"}) = read(seek(io, POINTDATA_OFFSET_SEEK), UInt32)

n_records(io::Stream{format"LAS"}) = read(seek(io, N_RECORDS_SEEK), UInt32)

#get indices of all fields that are not :x , :y, or :z
notcoord(ptype::Type{T}) where T <: AbstractPointDataRecordFormat = findall(Ref([:x, :y, :z]) .∌ fieldnames(ptype))

function calculate_offset(ptype::Type{T}, s::Symbol) where T <: AbstractPointDataRecordFormat
    if s == :coordinates
        offset = sum(sizeof.(ptype.types)[notcoord(ptype)])
    else
        offset = sum(sizeof.(ptype.types)[findall(!isequal(s), ptype.types)])
    end
    return offset
end

function get_transform(io::Stream{format"LAS"})
    seek(io, X_SCALE_SEEK)

    transform = AffineMap(Diagonal(SA_F64[read(io, Float64) 0 0; 0 read(io, Float64) 0; 0 0 read(io, Float64)]), SA_F64[read(io, Float64), read(io, Float64), read(io, Float64)])

    inverse_transform = (x->round.(Int32,x)) ∘ inv(transform)
    return transform, inverse_transform
end

function pointformat(io::Stream{format"LAS"})
    id = read(seek(io, FORMAT_SEEK), UInt8)
    if id == 0x00
        return PointDataRecordFormat0
    elseif id == 0x01
        return PointDataRecordFormat1
    elseif id == 0x02
        return PointDataRecordFormat2
    elseif id == 0x03
        return PointDataRecordFormat3
    else
        error("unsupported point format $(Int(id))")
    end
end

function load(f::File{format"LAS"}, field::Symbol)
    open(f) do s
        load(s, field)
    end
end

function load(s::Stream{format"LAS"}, field::Symbol)
    format = pointformat(s)
    records = n_records(s)
    point_seek = pointdata_start(s)
    if !(field in(fieldnames(format))) && field != :coordinates
        error("Field not part of LAS Point Data Records Format $format")
    elseif field == :coordinates
        coordinates = []
        transform,inverse = get_transform(s)
        seek(s, point_seek)
        while !eof(s)
            push!(coordinates, [read(s, Int32), read(s, Int32), read(s, Int32)])
            skip(s, calculate_offset(format, :coordinates))
        end
        coordinates_mapped = mappedarray(transform, inverse, coordinates)
        return coordinates_mapped
    else
        offset = calculate_offset(format, field)
        width = format.types[findfirst(isequal(field), fieldnames(format))]
        values = width[]
        seek(s, point_seek)
        while !eof(s)
            push!(values, read(s, width))
            skip(s, calculate_offset(format, field))
        end
        return values
    end
end

function __init__()
     add_format(format"LAS", "LASF", ".las", [:LASArrays])
end

end #end module


text:
using SurrogateModelOptim
using PlotlyJS
using Statistics

# Optimize the test function Rosenbrock
function rosenbrock_2D(x)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end
search_range=[(-5.0,5.0),(-5.0,5.0)]

result = smoptimize(rosenbrock_2D, search_range;
                    options=SurrogateModelOptim.Options(
                    SurrogateModelOptim.fast_opts(), #Run the suggested fast options to increase speed of surrogate optimisation
                    iterations=65,
                    create_final_surrogate=true, #Use the results from last iteration to
                                                 #re-create the surrogate before using it for plotting
                        ))
show(result)

function plot_fun_2D(fun,sr,title)
    N = 51
    x = range(sr[1][1], stop = sr[1][2], length = N)
    y = range(sr[2][1], stop = sr[2][2], length = N)
    grid(x,y) = [x,y]
    z = @. fun(grid(x,y'))

    trace = surface(x=x,y=y,z=z, colorscale="Viridis")
    layout = Layout(title=title)
    p = plot(trace,layout)
end

display(plot_fun_2D(rosenbrock_2D,search_range,"Original function"))
display(plot_fun_2D(x->median(result.sm_interpolant(x)),search_range,"Estimated function"))


text:
# Adjoint test for F and J
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall
using PyPlot
using opesciSLIM.TimeModeling

## Set up model structure
n = (160, 170)	# (x,y,z) or (x,z)
d = (15.,15.)
o = (0.,0.)

# Velocity [km/s]
v = ones(Float32,n) + 0.5f0
v[:,Int(round(end/3)):end] = 4f0
v0 = ones(Float32,n) + 0.5f0

# Slowness squared [s^2/km^2]
m = (1f0./v).^2
m0 = (1f0./v0).^2
dm = vec(m - m0)
rho =  v - .5f0
rho0 = v0 - .5f0

# Setup info and model structure
nsrc = 1
model = Model(n,d,o,m, rho)
model0 = Model(n,d,o,m0, rho0)

## Set up receiver geometry
nxrec = 160
xrec = linspace(100f0,1500f0,nxrec)
yrec = 0f0
zrec = linspace(30f0,30f0,nxrec)

# receiver sampling and recording time
timeR = 1500f0	# receiver recording time [ms]
dtR = 4.2f0	# receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = 1200f0
ysrc = 0f0
zsrc = 30f0

# source sampling and number of time steps
timeS = 1500f0
dtS = 4.2f0

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# Info structure
ntComp = get_computational_nt(srcGeometry,recGeometry,model)
info = Info(prod(n), nsrc, ntComp)

# setup wavelet
f0 = 0.01f0
wavelet = source(timeS,dtS,f0)
wave_rand = wavelet.*rand(Float32,size(wavelet))
###################################################################################################

# Modeling operators
F = joModeling(info,model,srcGeometry,recGeometry)
q = joData(srcGeometry,wavelet)
#
# Nonlinear modeling
d_hat = F*q

# Generate random noise data vector with size of d_hat
# in the range of F
qr = joData(srcGeometry,wave_rand)
d1 = F*qr

# Adjoint computation
q_hat = F'*d1

# Result F
println(abs(dot(d1,d_hat)))
println(abs(dot(q,q_hat)))
println("Residual: ", abs(dot(d1,d_hat) - dot(q,q_hat)))
println("Ratio: ", abs(dot(d1,d_hat)/dot(q,q_hat)))

# Linearized modeling
F0 = joModeling(info,model0,srcGeometry,recGeometry)
J = joJacobian(F0,q)

dD_hat = J*dm
dm_hat = J'*dD_hat

# Result J
println(dot(dD_hat,dD_hat))
println(dot(dm,dm_hat))
println("Residual: ", abs(dot(dD_hat,dD_hat) - dot(dm,dm_hat)))
println("Ratio: ", abs(dot(dD_hat,dD_hat)/dot(dm,dm_hat)))

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
module AllocationOpt

using CSV
using GraphQLClient

export network_state, optimize_indexer, read_filterlists, push_allocations!, create_rules!

include("exceptions.jl")
include("domainmodel.jl")
include("query.jl")
include("service.jl")
include("ActionQueue.jl")
include("CLI.jl")

"""
    function network_state(id, whitelist, blacklist, pinnedlist, frozenlist, indexer_service_network_url)

# Arguments
- `id::AbstractString`: The id of the indexer to optimise.
- `whitelist::Vector{AbstractString}`: Subgraph deployment IPFS hashes included in this list will be considered for, but not guaranteed allocation.
- `blacklist::Vector{AbstractString}`: Subgraph deployment IPFS hashes included in this list will not be considered, and will be suggested to close if there's an existing allocation.
- `pinnedlist::Vector{AbstractString}`: Subgraph deployment IPFS hashes included in this list will be guaranteed allocation. Currently unsupported.
- `frozenlist::Vector{AbstractString}`: Subgraph deployment IPFS hashes included in this list will not be considered during optimisation. Any allocations you have on these subgraphs deployments will remain.
- `indexer_service_network_url::AbstractString`: The URL that exposes the indexer service's network endpoint. Must begin with http. Example: http://localhost:7600/network.

"""
function network_state(
    id::AbstractString,
    whitelist::AbstractVector{T},
    blacklist::AbstractVector{T},
    pinnedlist::AbstractVector{T},
    frozenlist::AbstractVector{T},
    indexer_service_network_url::AbstractString,
) where {T<:AbstractString}
    if !isempty(pinnedlist)
        @warn "pinnedlist is not currently optimised for."
    end
    userlists = vcat(whitelist, blacklist, pinnedlist, frozenlist)
    if !verify_ipfshashes(userlists)
        throw(BadSubgraphIpfsHashError())
    end

    # Construct whitelist and blacklist
    query_ipfshash_in = ipfshash_in(whitelist, pinnedlist)
    query_ipfshash_not_in = ipfshash_not_in(blacklist, frozenlist)

    # Get client
    client = Client(indexer_service_network_url)

    # Pull data from mainnet subgraph
    repo = snapshot(client, query_ipfshash_in, query_ipfshash_not_in)

    # Handle frozenlist
    # Get indexer
    indexer, repo = detach_indexer(repo, id)

    # Reduce indexer stake by frozenlist
    fstake = frozen_stake(client, id, frozenlist)
    indexer = Indexer(indexer.id, indexer.stake - fstake, indexer.allocations)
    return repo, indexer
end

"""
    function optimize_indexer(indexer, repo, minimum_allocation_amount, maximum_new_allocations)

# Arguments
- `indexer::Indexer`: The indexer being optimised.
- `repo::Repository`: Contains the current network state.
- `minimum_allocation_amount::Real`: The minimum amount of GRT that you are willing to allocate to a subgraph.
- `maximum_new_allocations::Integer`: The maximum number of new allocations you would like the optimizer to open.
```
"""
function optimize_indexer(
    indexer::Indexer,
    repo::Repository,
    minimum_allocation_amount::Real,
    maximum_new_allocations::Integer,
)
    @warn "maximum_new_allocations is not currently optimised for."
    @warn "minimum_allocation_amount is not currently optimised for."

    # Optimise
    ω = optimize(indexer, repo)

    # Filter results with deployment IPFS hashes
    suggested_allocations = Dict(
        ipfshash(k) => v for (k, v) in zip(repo.subgraphs, ω) if v > 0.0
    )

    return suggested_allocations
end

"""
    function read_filterlists(filepath)

# Arguments

- `filepath::AbstractString`: A path to the CSV file that contains whitelist, blacklist, pinnedlist, frozenlist as columns.
"""
function read_filterlists(filepath::AbstractString)
    # Read file
    path = abspath(filepath)
    csv = CSV.File(path; header=1, types=String)

    # Filter out missings from csv
    listtypes = [:whitelist, :blacklist, :pinnedlist, :frozenlist]
    cols = map(x -> collect(skipmissing(csv[x])), listtypes)

    return cols
end

"""
    function push_allocations!(indexer_id, management_server_url, proposed_allocations, whitelist, blacklist, pinnedlist, frozenlist)

# Arguments

- `indexer_id::AbstractString`: Indexer id
- `management_server_url::T`: Indexer management server url, in format similar to http://localhost:18000
- `proposed_allocations::Dict{T,<:Real}`: The set of allocation to open returned by optimize_indexer
- `whitelist::AbstractVector{T}`: Unused. Here for completeness.
- `blacklist::AbstractVector{T}`: Unused. Here for completeness.
- `pinnedlist::AbstractVector{T}`: Unused. Here for completeness.
- `frozenlist::AbstractVector{T}`: Make sure to not close these allocations as they are frozen.
"""
function push_allocations!(
    indexer_id::AbstractString,
    management_server_url::T,
    indexer_service_network_url::T,
    proposed_allocations::Dict{T,<:Real},
    whitelist::AbstractVector{T},
    blacklist::AbstractVector{T},
    pinnedlist::AbstractVector{T},
    frozenlist::AbstractVector{T},
) where {T<:AbstractString}
    actions = []

    # Query existing allocations that are not frozen
    existing_allocations = query_indexer_allocations(
        Client(indexer_service_network_url), indexer_id
    )
    existing_allocs::Dict{String,String} = Dict(
        ipfshash.(existing_allocations) .=> id.(existing_allocations)
    )
    existing_ipfs::Vector{String} = ipfshash.(existing_allocations)
    proposed_ipfs = collect(keys(proposed_allocations))

    # Generate ActionQueue inputs
    reallocations, reallocate_ipfs = ActionQueue.reallocate_actions(
        proposed_ipfs, existing_ipfs, proposed_allocations, existing_allocs
    )
    open_allocations, open_ipfs = ActionQueue.allocate_actions(
        proposed_ipfs, reallocate_ipfs, proposed_allocations
    )
    close_allocations, close_ipfs = ActionQueue.unallocate_actions(
        existing_allocs, existing_ipfs, reallocate_ipfs, frozenlist
    )
    actions = vcat(reallocations, open_allocations, close_allocations)

    # Send ActionQueue inputs to indexer management server
    client = Client(management_server_url)
    response = mutate(client, "queueActions", Dict("actions" => actions))

    return response
end

function create_rules!(
    indexer_id::AbstractString,
    indexer_service_network_url::T,
    proposed_allocations::Dict{T,<:Real},
    whitelist::AbstractVector{T},
    blacklist::AbstractVector{T},
    pinnedlist::AbstractVector{T},
    frozenlist::AbstractVector{T},
) where {T<:AbstractString}
    actions = []

    # Query existing allocations that are not frozen
    existing_allocations = query_indexer_allocations(
        Client(indexer_service_network_url), indexer_id
    )
    existing_allocs = Dict(ipfshash.(existing_allocations) .=> id.(existing_allocations))
    existing_ipfs = ipfshash.(existing_allocations)
    proposed_ipfs = collect(keys(proposed_allocations))

    # Generate CLI commands

    reallocations, reallocate_ipfs = CLI.reallocate_actions(
        proposed_ipfs, existing_ipfs, proposed_allocations, existing_allocs
    )
    existing_ipfs ∩ proposed_ipfs
    open_allocations, open_ipfs = CLI.allocate_actions(
        proposed_ipfs, reallocate_ipfs, proposed_allocations
    )
    close_allocations, close_ipfs = CLI.unallocate_actions(
        existing_ipfs, reallocate_ipfs, frozenlist
    )
    actions = vcat(reallocations, open_allocations, close_allocations)
    return actions
end

end


text:
"""
uploads a texture by given file path
"""
function uploadTexture(path)
  img = Images.load(path)
  (imgwidth, imgheight) = size(img)
  imga = reinterpret.(vec(channelview(img)))

  texture = glGenTexture()
  glActiveTexture(GL_TEXTURE0)
  glBindTexture(GL_TEXTURE_2D, texture)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST) #GL_LINEAR
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST) #GL_LINEAR
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEAT
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEAT
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, imgwidth, imgheight, 0, GL_RGBA, GL_UNSIGNED_BYTE, imga)
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
using AlgebraicPetri
using Catlab.CategoricalAlgebra
using Catalyst
using OrdinaryDiffEq
using DiffEqFlux #=, Flux=#
using Plots
#import Catalyst: ReactionSystem

println("Starting")

#=SIR = LabelledPetriNet([:S, :I, :R],
  :inf => ((:S, :I)=>(:I, :I)),
  :rec => (:I=>:R),
  :id => (:S => :S),
  :id => (:I => :I),
  :id => (:R => :R)
)=#

get_infected_states(g::AbstractLabelledPetriNet) =
   [i for (i,s) in enumerate(g[:sname]) if occursin("I", string(s))]

get_susceptible_states(g::AbstractLabelledPetriNet) =
   [i for (i,s) in enumerate(g[:sname]) if occursin("S", string(s))]

get_recovered_states(g::AbstractLabelledPetriNet) =
   [i for (i,s) in enumerate(g[:sname]) if occursin("R", string(s)) || occursin("V", string(s))]

counter(a) = [count(==(i),a) for i in unique(a)]
function MakeReactionSystem(pn::AbstractPetriNet)
    @parameters t k[1:(nt(pn) + ns(pn))]
    @variables S[1:ns(pn)](t)

    rxs = map(1:nt(pn)) do t
      inpts = pn[incident(pn, t, :it),:is]
      otpts = pn[incident(pn, t, :ot),:os]
      in_count = collect(counter(inpts))
      ot_count = collect(counter(otpts))
      Reaction(k[t], [S[i] for i in unique(inpts)],
                     [S[o] for o in unique(otpts)],
                     in_count, ot_count)
    end

    ReactionSystem(rxs, t, S, k; name=:model)
end

#=SIR_rxn = MakeReactionSystem(SIR)
#AlgebraicPetri.Graph(SIR)
#Catalyst.Graph(SIR_rxn)

p_real = [.1/1000, .01]
tspan = (0.0,250.0)
u0 = [999.0, 1.0, 0.0]=#
#=sample_times = range(tspan[1],stop=tspan[2],length=100)

op = ODEProblem(SIR_rxn, u0, tspan, p_real)
sol_real = solve(op, Tsit5(), tstops=sample_times)

sample_vals = [sol_real.u[findfirst(sol_real.t .>= ts)][var] * (1+(0.1rand()-0.05))
                for var in 1:3, ts in sample_times]

plt = plot(sol_real, lw=2, label=reshape(map(string, SIR[:, :sname]), 1, ns(SIR)))
plot!(sample_times, sample_vals', seriestype=:scatter, label="")
display(plt)=#

"""
Generate infection data from a petri net model given:
  - rate parameters p
  - initial concentrations u0
  - timespan
  - number of samples
Returns sample infectious data, problem, and solution
"""
function generate_data(model::AbstractLabelledPetriNet, p, u0, tspan, num_samples)
    sample_times = range(tspan[1], stop=tspan[2], length=num_samples)

    prob = ODEProblem(MakeReactionSystem(model), u0, tspan, p)
    sol = solve(prob, Tsit5(), tstops=sample_times)

    inf_sample_vals = [sol.u[findfirst(sol.t .>= ts)][var] * (1+(0.1rand()-0.05))
                for var in get_infected_states(model), ts in sample_times]

    total_inf_samples = map(sum, eachcol(inf_sample_vals))

    susc_sample_vals = [sol.u[findfirst(sol.t .>= ts)][var] * (1+(0.1rand()-0.05))
                for var in get_susceptible_states(model), ts in sample_times]

    total_susc_samples = map(sum, eachcol(susc_sample_vals))

    rec_sample_vals = [sol.u[findfirst(sol.t .>= ts)][var] * (1+(0.1rand()-0.05))
                for var in get_recovered_states(model), ts in sample_times]

    total_rec_samples = map(sum, eachcol(rec_sample_vals))

    return hcat(total_inf_samples, total_rec_samples, total_susc_samples), sample_times, prob, sol
end

function optimise_p(model::AbstractLabelledPetriNet, op, p_init, tend, sample_data, sample_times)
    # Loss on all populations
    #=loss = function (p)
        sol = solve(remake(op,tspan=(0.0,tend),p=p), Tsit5(), tstops=sample_times)
        vals = hcat(map(ts -> sol.u[findfirst(sol.t .>= ts)], sample_times[1:findlast(sample_times .<= tend)])...)
        loss = sum(abs2, vals .- sample_vals[:,1:size(vals)[2]])
        return loss, sol
    end=#

    # Loss on just infected population (requires closer initial values)
    loss = function (p)
        sol = solve(remake(op,tspan=(0.0,tend),p=p), Tsit5(), tstops=sample_times)
        vals = hcat(map(ts -> sol.u[findfirst(sol.t .>= ts)], sample_times[1:findlast(sample_times .<= tend)])...)
        inf_vals = map(sum, collect(zip([vals[i,:] for i in get_infected_states(model)]...)))

        #loss = sum(abs2, vals[2,:] .- sample_vals[:,1:size(vals)[2]][2,:])
        inf_samples = sample_data[:,1]
        rec_samples = sample_data[:,2]
        susc_samples = sample_data[:,3]

        inf_loss = sum(abs2, inf_vals .- inf_samples[1:size(inf_vals)[1]])

        susc_vals = map(sum, collect(zip([vals[i,:] for i in get_susceptible_states(model)]...)))
        susc_loss = sum(abs2, susc_vals .- susc_samples[1:size(susc_vals)[1]])

        rec_vals = map(sum, collect(zip([vals[i,:] for i in get_recovered_states(model)]...)))
        rec_loss = sum(abs2, susc_vals .- rec_samples[1:size(rec_vals)[1]])
        return susc_loss + inf_loss, sol
    end

    callback = function (p, l, pred)
        display(l)
        plt = plot(pred, lw=2, linestyle=:dash, color=[:orange :green :blue])
        plot!(sample_times, sample_vals', seriestype=:scatter, label="")
        display(plt)
        return false
    end

    return DiffEqFlux.sciml_train(loss, p_init, #=ADAM(0.05), cb=callback,=# maxiters=500, abstol=1e-4, reltol=1e-4,
            lower_bounds=repeat([1e-6], length(p_init)), upper_bounds=ones(length(p_init)))
end

function full_train(model, u0, tspan, training_data, sample_times, param_guess)
    rxn = MakeReactionSystem(model)
    #p_estimate = zeros(numreactionparams(rxn))
    #p_estimate = repeat([1e-6], numreactionparams(rxn)-ns(model))
    p_estimate = param_guess
    prob = ODEProblem(rxn, u0, tspan, p_estimate)
    loss = 0
    sol_estimate = Nothing
    # TODO: algorithmically determine this range
    for i in 50:50:250
        res_ode = optimise_p(model, prob, p_estimate, i, training_data, sample_times)
        p_estimate = res_ode.minimizer
        sol_estimate = solve(remake(prob,tspan=tspan,p=p_estimate), Tsit5())

        susc_vals = map(sum, collect(zip([sol_estimate[i,:] for i in get_susceptible_states(model)]...)))
        rec_vals = map(sum, collect(zip([sol_estimate[i,:] for i in get_recovered_states(model)]...)))
        inf_vals = map(sum, collect(zip([sol_estimate[i,:] for i in get_infected_states(model)]...)))

        #=plt = plot(sample_times, training_data, seriestype=:scatter, label="")
        plot!(susc_vals, lw=2, label="S")
        plot!(rec_vals, lw=2, label="R")
        plot!(inf_vals, lw=2, label="I")
        display(plt)=#
        #=plt = plot(sample_times, training_data, seriestype=:scatter, label="")
        plot!(sol_estimate, lw=2, label=reshape(map(string, model[:, :sname]), 1, ns(model)))
        display(plt)=#

        println("Estimated params: $p_estimate")
        loss = res_ode.minimum
        println("Loss: $loss")
    end
    return  sol_estimate, loss
end

function stupid()
    return "dumb"
end

function full_train(model, u0, tspan, training_data, sample_times)
    rxn = MakeReactionSystem(model)
    #p_estimate = zeros(numreactionparams(rxn))
    p_estimate = repeat([1e-5], numreactionparams(rxn)-ns(model))
    #p_estimate = param_guess
    prob = ODEProblem(rxn, u0, tspan, p_estimate)
    loss = 0
    sol_estimate = Nothing
    # TODO: algorithmically determine this range
    for i in 50:50:250
        res_ode = optimise_p(model, prob, p_estimate, i, training_data, sample_times)
        p_estimate = res_ode.minimizer
        sol_estimate = solve(remake(prob,tspan=tspan,p=p_estimate), Tsit5())

        susc_vals = map(sum, collect(zip([sol_estimate[i,:] for i in get_susceptible_states(model)]...)))
        rec_vals = map(sum, collect(zip([sol_estimate[i,:] for i in get_recovered_states(model)]...)))
        inf_vals = map(sum, collect(zip([sol_estimate[i,:] for i in get_infected_states(model)]...)))

        #=plt = plot(sample_times, training_data, seriestype=:scatter, label="")
        plot!(susc_vals, lw=2, label="S")
        plot!(rec_vals, lw=2, label="R")
        plot!(inf_vals, lw=2, label="I")
        display(plt)=#
        #=plt = plot(sample_times, training_data, seriestype=:scatter, label="")
        plot!(sol_estimate, lw=2, label=reshape(map(string, model[:, :sname]), 1, ns(model)))
        display(plt)=#

        println("Estimated params: $p_estimate")
        loss = res_ode.minimum
        println("Loss: $loss")
    end
    return  sol_estimate, loss
end

function run_param_est()
    println("Generating data...")

    data, sample_times, op, sol_real = generate_data(SIR, p_real, u0, tspan, 100)
    plt = plot(sol_real, lw=2, label=reshape(map(string, SIR[:, :sname]), 1, ns(SIR)))
    plot!(sample_times, data, seriestype=:scatter, label="")
    display(plt)

    println("Running param estimation...")
    p_estimate = [0.0,0.0]
    for i in 50.0:50.0:250.0
        res_ode = optimise_p(SIR, op, p_estimate, i, data, sample_times)
        p_estimate = res_ode.minimizer
        sol_estimate = solve(remake(op,tspan=(0.0,250.0),p=p_estimate), Tsit5())

        plt = plot(sample_times, sample_vals', seriestype=:scatter, label="")
        plot!(sol_estimate, lw=2, label=reshape(map(string, SIR[:, :sname]), 1, ns(SIR)))
        display(plt)

        println("Estimated params: $p_estimate")
        loss = res_ode.minimum
        println("Loss: $loss")
    end

    sol_estimate = solve(remake(op,tspan=(0.0,250.0),p=p_estimate), Tsit5())
    return sol_estimate, p_estimate
end

#sol_estimate, p_estimate = run_param_est()

#=plot(sample_times, sample_vals', seriestype=:scatter, label="")
plot!(sol_estimate, lw=2)=#
#println("Final param estimate: $p_estimate")


text:
using MomentTensors, Test

@testset "Planes" begin
    @test all(isapprox.(MomentTensors.auxplane(30, 60, -45),
            (147.0, 52.0, -141.0), atol=1))
end


text:
struct Camera
    # Focal length.
    fx::Float64
    fy::Float64
    # Principal point.
    cx::Float64
    cy::Float64
    # Radial lens distortion.
    k1::Float64
    k2::Float64
    # Tangential lens distortion.
    p1::Float64
    p2::Float64
    # Intrinsic matrix.
    K::SMatrix{3, 3, Float64, 9}
    iK::SMatrix{3, 3, Float64, 9}
    # Image resolution.
    height::Int64
    width::Int64

    """
    Transformation from 0-th camera to i-th (this).
    """
    Ti0::SMatrix{4, 4, Float64, 16}
    """
    Transformation from i-th (this) camera to 0-th.
    """
    T0i::SMatrix{4, 4, Float64, 16}
end

function Camera(
    fx, fy, cx, cy, # Intrinsics.
    k1, k2, p1, p2, # Distortion coefficients.
    height, width; Ti0 = SMatrix{4, 4, Float64}(I),
)
    K = SMatrix{3, 3, Float64, 9}(
        fx, 0.0, 0.0,
        0.0, fy, 0.0,
        cx, cy, 1.0)
    Camera(
        fx, fy, cx, cy, k1, k2, p1, p2,
        K, inv(K), height, width, Ti0, inv(SE3, Ti0))
end

function Camera(;
    fx, fy, cx, cy, height, width,
    k1 = 0.0, k2 = 0.0, p1 = 0.0, p2 = 0.0,
    Ti0 = SMatrix{4, 4, Float64}(I),
)
    Camera(fx, fy, cx, cy, k1, k2, p1, p2, height, width; Ti0)
end

"""
Project point from 3D space onto the image plane.

# Arguments:
- `point`: Point in 3D space in `(x, y, z)` format.

# Returns:
Projected point in `(y, x)` format.
"""
function project(c::Camera, point)
    inv_z = 1.0 / point[3]
    Point2f(
        c.fy * point[2] * inv_z + c.cy,
        c.fx * point[1] * inv_z + c.cx)
end

"""
Project `point` onto image plane of the `Camera`,
accounting for the distortion parameters of the camera.

# Arguments:
- `point`: 3D point to project in `(x, y, z)` format.

# Returns:
2D floating point coordinates in `(y, x)` format.
"""
function project_undistort(c::Camera, point)
    normalized = point[2:-1:1] ./ point[3]
    undistort_pdn_point(c, normalized)
end

"""
Check if `point` is in the image bounds of the `Camera`.

# Arguments:
- `point`: Point to check. In `(y, x)` format.
"""
function in_image(c::Camera, point)
    1 ≤ point[1] ≤ c.height && 1 ≤ point[2] ≤ c.width
end

"""
# Arguments:
- `point::SVector{2}`: Point to undistort. In `(y, x)` format.
"""
function undistort_point(c::Camera, point::Point2)
    normalized = Point2f(
        (point[1] - c.cy) / c.fy,
        (point[2] - c.cx) / c.fx)
    undistort_pdn_point(c, normalized)
end

"""
Undistort point.

# Arguments:
- `point::SVector{2}`: Predivided by `K` & normalized point in `(y, x)` format.
"""
function undistort_pdn_point(c::Camera, point)
    sqrd_normalized = point .^ 2
    # Square radius from center.
    sqrd_radius = sqrd_normalized |> sum
    # Radial distortion factor.
    rd = 1.0 + c.k1 * sqrd_radius + c.k2 * sqrd_radius ^ 2
    # Tangential distortion component.
    p = point |> prod
    dtx = 2 * c.p1 * p + c.p2 * (sqrd_radius + 2 * sqrd_normalized[1])
    dty = c.p1 * (sqrd_radius + 2 * sqrd_normalized[2]) + 2 * c.p2 * p
    # Lens distortion coordinates.
    distorted = rd .* point .+ (dty, dtx)
    # Final projection (assume skew is always `0`).
    Point2f(distorted .* (c.fy, c.fx) .+ (c.cy, c.cx))
end

"""
Transform point from 2D to 3D by dividing by `K`.

# Arguments:

- `point::Point2`: Point to backproject in `(y, x)` format.

# Returns:

Backprojected `Point3f` in `(x, y, z = 1.0)` format.
"""
function backproject(c::Camera, point::Point2)
    Point3f((point[2] - c.cx) / c.fx, (point[1] - c.cy) / c.fy, 1.0)
end


text:
module test_bukdu_system

module NAStuff

using Bukdu
struct NA <: ApplicationController; conn::Conn end
index(c::NA) = nothing
hello(c::NA) = TEST_INTERNAL_ERROR

end # module NAStuff

using Test
using Bukdu
using .Bukdu.System

function index
end

routes() do
    get("/na", NAStuff.NA, index)
    post("/hello", NAStuff.NA, NAStuff.hello)
end

result = Router.call(get, "/na")
@test occursin("Bukdu.System.NotApplicableError", result.got)
@test result.resp.status == 500
@test result.route.action === System.not_applicable

Plug.Loggers.config[:error_stackframes_range] = 1:2
result = Router.call(post, "/hello")
@test result.route.action === System.internal_error

Routing.empty!()

@test Plug.Loggers._regularize_text("가1", 1) == "가"

struct Controller <: ApplicationController
    conn::Conn
end
function index(::Controller)
end

struct VeryLongNamedController <: ApplicationController
    conn::Conn
end
function index(::VeryLongNamedController)
end

routes() do
    get("/just", Controller, index)
    get("/long", VeryLongNamedController, index)
end
Router.call(get, "/just")
Router.call(get, "/long")

Routing.empty!()

end # module test_bukdu_system