text:
# Julia wrapper for header: /usr/include/scip/pub_var.h
# Automatically generated using Clang.jl wrap_c

function SCIPvarGetNLocksDown(var)
    ccall((:SCIPvarGetNLocksDown, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNLocksUp(var)
    ccall((:SCIPvarGetNLocksUp, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNLocksUpType(var, locktype)
    ccall((:SCIPvarGetNLocksUpType, libscip), Cint, (Ptr{SCIP_VAR}, SCIP_LOCKTYPE), var, locktype)
end

function SCIPvarGetNLocksDownType(var, locktype)
    ccall((:SCIPvarGetNLocksDownType, libscip), Cint, (Ptr{SCIP_VAR}, SCIP_LOCKTYPE), var, locktype)
end

function SCIPvarMayRoundDown(var)
    ccall((:SCIPvarMayRoundDown, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarMayRoundUp(var)
    ccall((:SCIPvarMayRoundUp, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarCompareActiveAndNegated(var1, var2)
    ccall((:SCIPvarCompareActiveAndNegated, libscip), Cint, (Ptr{SCIP_VAR}, Ptr{SCIP_VAR}), var1, var2)
end

function SCIPvarCompActiveAndNegated(elem1, elem2)
    ccall((:SCIPvarCompActiveAndNegated, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPvarCompare(var1, var2)
    ccall((:SCIPvarCompare, libscip), Cint, (Ptr{SCIP_VAR}, Ptr{SCIP_VAR}), var1, var2)
end

function SCIPvarComp(elem1, elem2)
    ccall((:SCIPvarComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPvarCompObj(elem1, elem2)
    ccall((:SCIPvarCompObj, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPvarGetHashkey(userptr, elem)
    ccall((:SCIPvarGetHashkey, libscip), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), userptr, elem)
end

function SCIPvarIsHashkeyEq(userptr, key1, key2)
    ccall((:SCIPvarIsHashkeyEq, libscip), UInt32, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), userptr, key1, key2)
end

function SCIPvarGetHashkeyVal(userptr, key)
    ccall((:SCIPvarGetHashkeyVal, libscip), UInt64, (Ptr{Cvoid}, Ptr{Cvoid}), userptr, key)
end

function SCIPvarsGetProbvar(vars, nvars)
    ccall((:SCIPvarsGetProbvar, libscip), Cvoid, (Ptr{Ptr{SCIP_VAR}}, Cint), vars, nvars)
end

function SCIPvarGetProbvar(var)
    ccall((:SCIPvarGetProbvar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarsGetProbvarBinary(vars, negatedarr, nvars)
    ccall((:SCIPvarsGetProbvarBinary, libscip), SCIP_RETCODE, (Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Ptr{UInt32}}, Cint), vars, negatedarr, nvars)
end

function SCIPvarGetProbvarBinary(var, negated)
    ccall((:SCIPvarGetProbvarBinary, libscip), SCIP_RETCODE, (Ptr{Ptr{SCIP_VAR}}, Ptr{UInt32}), var, negated)
end

function SCIPvarGetProbvarBound(var, bound, boundtype)
    ccall((:SCIPvarGetProbvarBound, libscip), SCIP_RETCODE, (Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{SCIP_BOUNDTYPE}), var, bound, boundtype)
end

function SCIPvarGetProbvarHole(var, left, right)
    ccall((:SCIPvarGetProbvarHole, libscip), SCIP_RETCODE, (Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{Cdouble}), var, left, right)
end

function SCIPvarGetOrigvarSum(var, scalar, constant)
    ccall((:SCIPvarGetOrigvarSum, libscip), SCIP_RETCODE, (Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{Cdouble}), var, scalar, constant)
end

function SCIPvarIsTransformedOrigvar(var)
    ccall((:SCIPvarIsTransformedOrigvar, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNBranchings(var, dir)
    ccall((:SCIPvarGetNBranchings, libscip), Clonglong, (Ptr{SCIP_VAR}, SCIP_BRANCHDIR), var, dir)
end

function SCIPvarGetNBranchingsCurrentRun(var, dir)
    ccall((:SCIPvarGetNBranchingsCurrentRun, libscip), Clonglong, (Ptr{SCIP_VAR}, SCIP_BRANCHDIR), var, dir)
end

function SCIPvarGetInferenceSum(var, dir)
    ccall((:SCIPvarGetInferenceSum, libscip), Cdouble, (Ptr{SCIP_VAR}, SCIP_BRANCHDIR), var, dir)
end

function SCIPvarGetInferenceSumCurrentRun(var, dir)
    ccall((:SCIPvarGetInferenceSumCurrentRun, libscip), Cdouble, (Ptr{SCIP_VAR}, SCIP_BRANCHDIR), var, dir)
end

function SCIPvarGetCutoffSum(var, dir)
    ccall((:SCIPvarGetCutoffSum, libscip), Cdouble, (Ptr{SCIP_VAR}, SCIP_BRANCHDIR), var, dir)
end

function SCIPvarGetCutoffSumCurrentRun(var, dir)
    ccall((:SCIPvarGetCutoffSumCurrentRun, libscip), Cdouble, (Ptr{SCIP_VAR}, SCIP_BRANCHDIR), var, dir)
end

function SCIPvarGetAvgBranchdepth(var, dir)
    ccall((:SCIPvarGetAvgBranchdepth, libscip), Cdouble, (Ptr{SCIP_VAR}, SCIP_BRANCHDIR), var, dir)
end

function SCIPvarGetAvgBranchdepthCurrentRun(var, dir)
    ccall((:SCIPvarGetAvgBranchdepthCurrentRun, libscip), Cdouble, (Ptr{SCIP_VAR}, SCIP_BRANCHDIR), var, dir)
end

function SCIPvarHasImplic(var, varfixing, implvar, impltype)
    ccall((:SCIPvarHasImplic, libscip), UInt32, (Ptr{SCIP_VAR}, UInt32, Ptr{SCIP_VAR}, SCIP_BOUNDTYPE), var, varfixing, implvar, impltype)
end

function SCIPvarHasBinaryImplic(var, varfixing, implvar, implvarfixing)
    ccall((:SCIPvarHasBinaryImplic, libscip), UInt32, (Ptr{SCIP_VAR}, UInt32, Ptr{SCIP_VAR}, UInt32), var, varfixing, implvar, implvarfixing)
end

function SCIPvarsHaveCommonClique(var1, value1, var2, value2, regardimplics)
    ccall((:SCIPvarsHaveCommonClique, libscip), UInt32, (Ptr{SCIP_VAR}, UInt32, Ptr{SCIP_VAR}, UInt32, UInt32), var1, value1, var2, value2, regardimplics)
end

function SCIPvarGetAggregatedObj(var, aggrobj)
    ccall((:SCIPvarGetAggregatedObj, libscip), SCIP_RETCODE, (Ptr{SCIP_VAR}, Ptr{Cdouble}), var, aggrobj)
end

function SCIPvarSetInitial(var, initial)
    ccall((:SCIPvarSetInitial, libscip), SCIP_RETCODE, (Ptr{SCIP_VAR}, UInt32), var, initial)
end

function SCIPvarSetRemovable(var, removable)
    ccall((:SCIPvarSetRemovable, libscip), SCIP_RETCODE, (Ptr{SCIP_VAR}, UInt32), var, removable)
end

function SCIPvarGetName(var)
    ccall((:SCIPvarGetName, libscip), Cstring, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNUses(var)
    ccall((:SCIPvarGetNUses, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetData(var)
    ccall((:SCIPvarGetData, libscip), Ptr{SCIP_VARDATA}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarSetData(var, vardata)
    ccall((:SCIPvarSetData, libscip), Cvoid, (Ptr{SCIP_VAR}, Ptr{SCIP_VARDATA}), var, vardata)
end

function SCIPvarSetDelorigData(var, vardelorig)
    ccall((:SCIPvarSetDelorigData, libscip), Cvoid, (Ptr{SCIP_VAR}, Ptr{Cvoid}), var, vardelorig)
end

function SCIPvarSetTransData(var, vartrans)
    ccall((:SCIPvarSetTransData, libscip), Cvoid, (Ptr{SCIP_VAR}, Ptr{Cvoid}), var, vartrans)
end

function SCIPvarSetDeltransData(var, vardeltrans)
    ccall((:SCIPvarSetDeltransData, libscip), Cvoid, (Ptr{SCIP_VAR}, Ptr{Cvoid}), var, vardeltrans)
end

function SCIPvarSetCopyData(var, varcopy)
    ccall((:SCIPvarSetCopyData, libscip), Cvoid, (Ptr{SCIP_VAR}, Ptr{Cvoid}), var, varcopy)
end

function SCIPvarGetStatus(var)
    ccall((:SCIPvarGetStatus, libscip), SCIP_VARSTATUS, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsOriginal(var)
    ccall((:SCIPvarIsOriginal, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsTransformed(var)
    ccall((:SCIPvarIsTransformed, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsNegated(var)
    ccall((:SCIPvarIsNegated, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetType(var)
    ccall((:SCIPvarGetType, libscip), SCIP_VARTYPE, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsBinary(var)
    ccall((:SCIPvarIsBinary, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsIntegral(var)
    ccall((:SCIPvarIsIntegral, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsInitial(var)
    ccall((:SCIPvarIsInitial, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsRemovable(var)
    ccall((:SCIPvarIsRemovable, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsDeleted(var)
    ccall((:SCIPvarIsDeleted, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarMarkDeletable(var)
    ccall((:SCIPvarMarkDeletable, libscip), Cvoid, (Ptr{SCIP_VAR},), var)
end

function SCIPvarMarkNotDeletable(var)
    ccall((:SCIPvarMarkNotDeletable, libscip), Cvoid, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsDeletable(var)
    ccall((:SCIPvarIsDeletable, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarMarkDeleteGlobalStructures(var)
    ccall((:SCIPvarMarkDeleteGlobalStructures, libscip), Cvoid, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsActive(var)
    ccall((:SCIPvarIsActive, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetIndex(var)
    ccall((:SCIPvarGetIndex, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetProbindex(var)
    ccall((:SCIPvarGetProbindex, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetTransVar(var)
    ccall((:SCIPvarGetTransVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetCol(var)
    ccall((:SCIPvarGetCol, libscip), Ptr{SCIP_COL}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarIsInLP(var)
    ccall((:SCIPvarIsInLP, libscip), UInt32, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetAggrVar(var)
    ccall((:SCIPvarGetAggrVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetAggrScalar(var)
    ccall((:SCIPvarGetAggrScalar, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetAggrConstant(var)
    ccall((:SCIPvarGetAggrConstant, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetMultaggrNVars(var)
    ccall((:SCIPvarGetMultaggrNVars, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetMultaggrVars(var)
    ccall((:SCIPvarGetMultaggrVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetMultaggrScalars(var)
    ccall((:SCIPvarGetMultaggrScalars, libscip), Ptr{Cdouble}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetMultaggrConstant(var)
    ccall((:SCIPvarGetMultaggrConstant, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNegatedVar(var)
    ccall((:SCIPvarGetNegatedVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNegationVar(var)
    ccall((:SCIPvarGetNegationVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNegationConstant(var)
    ccall((:SCIPvarGetNegationConstant, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetObj(var)
    ccall((:SCIPvarGetObj, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetUnchangedObj(var)
    ccall((:SCIPvarGetUnchangedObj, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetLbOriginal(var)
    ccall((:SCIPvarGetLbOriginal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetUbOriginal(var)
    ccall((:SCIPvarGetUbOriginal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetHolelistOriginal(var)
    ccall((:SCIPvarGetHolelistOriginal, libscip), Ptr{SCIP_HOLELIST}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetLbGlobal(var)
    ccall((:SCIPvarGetLbGlobal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetUbGlobal(var)
    ccall((:SCIPvarGetUbGlobal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetHolelistGlobal(var)
    ccall((:SCIPvarGetHolelistGlobal, libscip), Ptr{SCIP_HOLELIST}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBestBoundGlobal(var)
    ccall((:SCIPvarGetBestBoundGlobal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetWorstBoundGlobal(var)
    ccall((:SCIPvarGetWorstBoundGlobal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetLbLocal(var)
    ccall((:SCIPvarGetLbLocal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetUbLocal(var)
    ccall((:SCIPvarGetUbLocal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetHolelistLocal(var)
    ccall((:SCIPvarGetHolelistLocal, libscip), Ptr{SCIP_HOLELIST}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBestBoundLocal(var)
    ccall((:SCIPvarGetBestBoundLocal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetWorstBoundLocal(var)
    ccall((:SCIPvarGetWorstBoundLocal, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBestBoundType(var)
    ccall((:SCIPvarGetBestBoundType, libscip), SCIP_BOUNDTYPE, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetWorstBoundType(var)
    ccall((:SCIPvarGetWorstBoundType, libscip), SCIP_BOUNDTYPE, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetLbLazy(var)
    ccall((:SCIPvarGetLbLazy, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetUbLazy(var)
    ccall((:SCIPvarGetUbLazy, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBranchFactor(var)
    ccall((:SCIPvarGetBranchFactor, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBranchPriority(var)
    ccall((:SCIPvarGetBranchPriority, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBranchDirection(var)
    ccall((:SCIPvarGetBranchDirection, libscip), SCIP_BRANCHDIR, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNVlbs(var)
    ccall((:SCIPvarGetNVlbs, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetVlbVars(var)
    ccall((:SCIPvarGetVlbVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetVlbCoefs(var)
    ccall((:SCIPvarGetVlbCoefs, libscip), Ptr{Cdouble}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetVlbConstants(var)
    ccall((:SCIPvarGetVlbConstants, libscip), Ptr{Cdouble}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNVubs(var)
    ccall((:SCIPvarGetNVubs, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetVubVars(var)
    ccall((:SCIPvarGetVubVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetVubCoefs(var)
    ccall((:SCIPvarGetVubCoefs, libscip), Ptr{Cdouble}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetVubConstants(var)
    ccall((:SCIPvarGetVubConstants, libscip), Ptr{Cdouble}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNImpls(var, varfixing)
    ccall((:SCIPvarGetNImpls, libscip), Cint, (Ptr{SCIP_VAR}, UInt32), var, varfixing)
end

function SCIPvarGetImplVars(var, varfixing)
    ccall((:SCIPvarGetImplVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_VAR}, UInt32), var, varfixing)
end

function SCIPvarGetImplTypes(var, varfixing)
    ccall((:SCIPvarGetImplTypes, libscip), Ptr{SCIP_BOUNDTYPE}, (Ptr{SCIP_VAR}, UInt32), var, varfixing)
end

function SCIPvarGetImplBounds(var, varfixing)
    ccall((:SCIPvarGetImplBounds, libscip), Ptr{Cdouble}, (Ptr{SCIP_VAR}, UInt32), var, varfixing)
end

function SCIPvarGetImplIds(var, varfixing)
    ccall((:SCIPvarGetImplIds, libscip), Ptr{Cint}, (Ptr{SCIP_VAR}, UInt32), var, varfixing)
end

function SCIPvarGetNCliques(var, varfixing)
    ccall((:SCIPvarGetNCliques, libscip), Cint, (Ptr{SCIP_VAR}, UInt32), var, varfixing)
end

function SCIPvarGetCliques(var, varfixing)
    ccall((:SCIPvarGetCliques, libscip), Ptr{Ptr{SCIP_CLIQUE}}, (Ptr{SCIP_VAR}, UInt32), var, varfixing)
end

function SCIPvarGetLPSol(var)
    ccall((:SCIPvarGetLPSol, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNLPSol(var)
    ccall((:SCIPvarGetNLPSol, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBdchgInfoLb(var, pos)
    ccall((:SCIPvarGetBdchgInfoLb, libscip), Ptr{SCIP_BDCHGINFO}, (Ptr{SCIP_VAR}, Cint), var, pos)
end

function SCIPvarGetNBdchgInfosLb(var)
    ccall((:SCIPvarGetNBdchgInfosLb, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBdchgInfoUb(var, pos)
    ccall((:SCIPvarGetBdchgInfoUb, libscip), Ptr{SCIP_BDCHGINFO}, (Ptr{SCIP_VAR}, Cint), var, pos)
end

function SCIPvarGetNBdchgInfosUb(var)
    ccall((:SCIPvarGetNBdchgInfosUb, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetValuehistory(var)
    ccall((:SCIPvarGetValuehistory, libscip), Ptr{SCIP_VALUEHISTORY}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetLPSol_rec(var)
    ccall((:SCIPvarGetLPSol_rec, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetNLPSol_rec(var)
    ccall((:SCIPvarGetNLPSol_rec, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetPseudoSol(var)
    ccall((:SCIPvarGetPseudoSol, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetSol(var, getlpval)
    ccall((:SCIPvarGetSol, libscip), Cdouble, (Ptr{SCIP_VAR}, UInt32), var, getlpval)
end

function SCIPvarGetRootSol(var)
    ccall((:SCIPvarGetRootSol, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBestRootSol(var)
    ccall((:SCIPvarGetBestRootSol, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBestRootRedcost(var)
    ccall((:SCIPvarGetBestRootRedcost, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetBestRootLPObjval(var)
    ccall((:SCIPvarGetBestRootLPObjval, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarSetBestRootSol(var, rootsol, rootredcost, rootlpobjval)
    ccall((:SCIPvarSetBestRootSol, libscip), Cvoid, (Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble), var, rootsol, rootredcost, rootlpobjval)
end

function SCIPvarGetAvgSol(var)
    ccall((:SCIPvarGetAvgSol, libscip), Cdouble, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetLbchgInfo(var, bdchgidx, after)
    ccall((:SCIPvarGetLbchgInfo, libscip), Ptr{SCIP_BDCHGINFO}, (Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, UInt32), var, bdchgidx, after)
end

function SCIPvarGetUbchgInfo(var, bdchgidx, after)
    ccall((:SCIPvarGetUbchgInfo, libscip), Ptr{SCIP_BDCHGINFO}, (Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, UInt32), var, bdchgidx, after)
end

function SCIPvarGetBdchgInfo(var, boundtype, bdchgidx, after)
    ccall((:SCIPvarGetBdchgInfo, libscip), Ptr{SCIP_BDCHGINFO}, (Ptr{SCIP_VAR}, SCIP_BOUNDTYPE, Ptr{SCIP_BDCHGIDX}, UInt32), var, boundtype, bdchgidx, after)
end

function SCIPvarGetLbAtIndex(var, bdchgidx, after)
    ccall((:SCIPvarGetLbAtIndex, libscip), Cdouble, (Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, UInt32), var, bdchgidx, after)
end

function SCIPvarGetUbAtIndex(var, bdchgidx, after)
    ccall((:SCIPvarGetUbAtIndex, libscip), Cdouble, (Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, UInt32), var, bdchgidx, after)
end

function SCIPvarGetBdAtIndex(var, boundtype, bdchgidx, after)
    ccall((:SCIPvarGetBdAtIndex, libscip), Cdouble, (Ptr{SCIP_VAR}, SCIP_BOUNDTYPE, Ptr{SCIP_BDCHGIDX}, UInt32), var, boundtype, bdchgidx, after)
end

function SCIPvarWasFixedAtIndex(var, bdchgidx, after)
    ccall((:SCIPvarWasFixedAtIndex, libscip), UInt32, (Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, UInt32), var, bdchgidx, after)
end

function SCIPvarGetLastBdchgIndex(var)
    ccall((:SCIPvarGetLastBdchgIndex, libscip), Ptr{SCIP_BDCHGIDX}, (Ptr{SCIP_VAR},), var)
end

function SCIPvarGetLastBdchgDepth(var)
    ccall((:SCIPvarGetLastBdchgDepth, libscip), Cint, (Ptr{SCIP_VAR},), var)
end

function SCIPvarWasFixedEarlier(var1, var2)
    ccall((:SCIPvarWasFixedEarlier, libscip), UInt32, (Ptr{SCIP_VAR}, Ptr{SCIP_VAR}), var1, var2)
end

function SCIPbdchgidxIsEarlier(bdchgidx1, bdchgidx2)
    ccall((:SCIPbdchgidxIsEarlier, libscip), UInt32, (Ptr{SCIP_BDCHGIDX}, Ptr{SCIP_BDCHGIDX}), bdchgidx1, bdchgidx2)
end

function SCIPbdchgidxIsEarlierNonNull(bdchgidx1, bdchgidx2)
    ccall((:SCIPbdchgidxIsEarlierNonNull, libscip), UInt32, (Ptr{SCIP_BDCHGIDX}, Ptr{SCIP_BDCHGIDX}), bdchgidx1, bdchgidx2)
end

function SCIPbdchginfoGetOldbound(bdchginfo)
    ccall((:SCIPbdchginfoGetOldbound, libscip), Cdouble, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetNewbound(bdchginfo)
    ccall((:SCIPbdchginfoGetNewbound, libscip), Cdouble, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetVar(bdchginfo)
    ccall((:SCIPbdchginfoGetVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetChgtype(bdchginfo)
    ccall((:SCIPbdchginfoGetChgtype, libscip), SCIP_BOUNDCHGTYPE, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetBoundtype(bdchginfo)
    ccall((:SCIPbdchginfoGetBoundtype, libscip), SCIP_BOUNDTYPE, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetDepth(bdchginfo)
    ccall((:SCIPbdchginfoGetDepth, libscip), Cint, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetPos(bdchginfo)
    ccall((:SCIPbdchginfoGetPos, libscip), Cint, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetIdx(bdchginfo)
    ccall((:SCIPbdchginfoGetIdx, libscip), Ptr{SCIP_BDCHGIDX}, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetInferVar(bdchginfo)
    ccall((:SCIPbdchginfoGetInferVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetInferCons(bdchginfo)
    ccall((:SCIPbdchginfoGetInferCons, libscip), Ptr{SCIP_CONS}, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetInferProp(bdchginfo)
    ccall((:SCIPbdchginfoGetInferProp, libscip), Ptr{SCIP_PROP}, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetInferInfo(bdchginfo)
    ccall((:SCIPbdchginfoGetInferInfo, libscip), Cint, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoGetInferBoundtype(bdchginfo)
    ccall((:SCIPbdchginfoGetInferBoundtype, libscip), SCIP_BOUNDTYPE, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoIsRedundant(bdchginfo)
    ccall((:SCIPbdchginfoIsRedundant, libscip), UInt32, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoHasInferenceReason(bdchginfo)
    ccall((:SCIPbdchginfoHasInferenceReason, libscip), UInt32, (Ptr{SCIP_BDCHGINFO},), bdchginfo)
end

function SCIPbdchginfoIsTighter(bdchginfo1, bdchginfo2)
    ccall((:SCIPbdchginfoIsTighter, libscip), UInt32, (Ptr{SCIP_BDCHGINFO}, Ptr{SCIP_BDCHGINFO}), bdchginfo1, bdchginfo2)
end

function SCIPboundchgGetNewbound(boundchg)
    ccall((:SCIPboundchgGetNewbound, libscip), Cdouble, (Ptr{SCIP_BOUNDCHG},), boundchg)
end

function SCIPboundchgGetVar(boundchg)
    ccall((:SCIPboundchgGetVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_BOUNDCHG},), boundchg)
end

function SCIPboundchgGetBoundchgtype(boundchg)
    ccall((:SCIPboundchgGetBoundchgtype, libscip), SCIP_BOUNDCHGTYPE, (Ptr{SCIP_BOUNDCHG},), boundchg)
end

function SCIPboundchgGetBoundtype(boundchg)
    ccall((:SCIPboundchgGetBoundtype, libscip), SCIP_BOUNDTYPE, (Ptr{SCIP_BOUNDCHG},), boundchg)
end

function SCIPboundchgIsRedundant(boundchg)
    ccall((:SCIPboundchgIsRedundant, libscip), UInt32, (Ptr{SCIP_BOUNDCHG},), boundchg)
end

function SCIPdomchgGetNBoundchgs(domchg)
    ccall((:SCIPdomchgGetNBoundchgs, libscip), Cint, (Ptr{SCIP_DOMCHG},), domchg)
end

function SCIPdomchgGetBoundchg(domchg, pos)
    ccall((:SCIPdomchgGetBoundchg, libscip), Ptr{SCIP_BOUNDCHG}, (Ptr{SCIP_DOMCHG}, Cint), domchg, pos)
end

function SCIPholelistGetLeft(holelist)
    ccall((:SCIPholelistGetLeft, libscip), Cdouble, (Ptr{SCIP_HOLELIST},), holelist)
end

function SCIPholelistGetRight(holelist)
    ccall((:SCIPholelistGetRight, libscip), Cdouble, (Ptr{SCIP_HOLELIST},), holelist)
end

function SCIPholelistGetNext(holelist)
    ccall((:SCIPholelistGetNext, libscip), Ptr{SCIP_HOLELIST}, (Ptr{SCIP_HOLELIST},), holelist)
end


text:
using KitBase, Plots, OffsetArrays
using ProgressMeter: @showprogress
using Base.Threads: @threads

set = Setup(space = "1d1f1v", boundary = "fix", maxTime = 0.1)
ps = PSpace1D(0, 1, 100, 1)
vs = VSpace1D(-5, 5, 28)
gas = Gas(Kn = 5e-4, K = 0.0, γ = 3.0)
ib = begin
    _fw = function (x)
        # wave
        #ρ = 1 + 0.1 * sin(2π * x)
        #U = 1.0
        #T = 1 / ρ

        # sod
        ρ, U, T = begin
            if x < 0.5
                1.0, 0.0, 1 / 0.5
            else
                0.125, 0.0, 1 / 0.625
            end
        end
        prim = [ρ, U, 1 / T]
        return prim_conserve(prim, gas.γ)
    end
    _ff = function (x)
        w = _fw(x)
        prim = conserve_prim(w, gas.γ)
        return maxwellian(vs.u, prim)
    end
    _bc = function (x)
        w = _fw(x)
        return conserve_prim(w, gas.γ)
    end
    IB1F(_fw, _ff, _bc)
end
ks = SolverSet(set, ps, vs, gas, ib)
#=
ctr, face = init_fvm(ks)
plot(ks, ctr)

t = 0.0
dt = timestep(ks, ctr, t)
nt = ks.set.maxTime ÷ dt |> Int
@showprogress for iter = 1:nt
    evolve!(ks, ctr, face, dt)
    update!(ks, ctr, face, dt, zeros(3))
end

plot!(ks, ctr)
=#

# basis
#m = zeros(5, length(vs.u))
m = zeros(7, length(vs.u))
for j in axes(m, 2)
    for i in axes(m, 1)
        m[i, j] = vs.u[j]^(i - 1)
    end
end
α = zeros(size(m, 1)) # minimizer

###
# one-cell test
###
#f0 = maxwellian(vs.u, [1.0, 0.0, 1.0])
#=f0 = 0.4 .* maxwellian(vs.u, [1.0, 1.0, 1.0]) .+
    0.6 .* maxwellian(vs.u, [1.0, -1.0, 1.0])

u0 = [sum(vs.weights .* m[i, :] .* f0) for i in axes(m, 1)]
α = zeros(size(m, 1))
res = KitBase.optimize_closure(α, m, vs.weights, u0, exp)
f = exp.(res.minimizer' * m)'

plot(vs.u, f0)
plot!(vs.u, f)=#

ctr = OffsetArray{ControlVolume1F}(undef, eachindex(ks.ps.x))
for i in eachindex(ctr)
    cons = ks.ib.fw(ks.ps.x[i])
    prim = conserve_prim(cons, ks.gas.γ)
    f = ks.ib.ff(ks.ps.x[i])

    w = zeros(size(m, 1))
    for j in axes(w, 1)
        w[j] = sum(ks.vs.weights .* m[j, :] .* f)
    end

    ctr[i] = ControlVolume(w, prim, f, 1)
end

face = Array{Interface1F}(undef, ks.ps.nx + 1)
for i = 1:ks.ps.nx+1
    fw = zeros(size(m, 1))
    ff = zero(ks.vs.u)
    face[i] = Interface(fw, ff, 1)
end

function rc!(ks, ctr, p)
    m, α = p

    @inbounds @threads for i = 0:ks.ps.nx+1
        res = KitBase.optimize_closure(α, m, ks.vs.weights, ctr[i].w, exp)
        ctr[i].f .= exp.(res.minimizer' * m)'
    end
end

function ev!(ks, ctr, face, p)
    m, dt = p

    @inbounds @threads for i = 1:ks.ps.nx+1
        flux_kfvs!(face[i].ff, ctr[i-1].f, ctr[i].f, ks.vs.u, dt)
        for j in axes(face[i].fw, 1)
            face[i].fw[j] = sum(ks.vs.weights .* m[j, :] .* face[i].ff)
        end
    end
end

function up!(ks, ctr, face, p)
    m, dt = p

    @inbounds @threads for i = 1:ks.ps.nx
        M = maxwellian(ks.vs.u, ctr[i].prim)
        τ = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
        q = @. (M - ctr[i].f) / τ

        Q = zeros(size(m, 1))
        for j = 4:size(m, 1)
            @inbounds Q[j] = sum(ks.vs.weights .* m[j, :] .* q)
        end

        @. ctr[i].w += (face[i].fw - face[i+1].fw) / ks.ps.dx[i] + Q * dt
        ctr[i].prim .= conserve_prim([ctr[i].w[1:2]; ctr[i].w[3] / 2], ks.gas.γ)
    end
end

t = 0.0
dt = timestep(ks, ctr, t)
nt = ks.set.maxTime ÷ dt |> Int

@showprogress for iter = 1:nt
    rc!(ks, ctr, (m, α))
    ev!(ks, ctr, face, (m, dt))
    up!(ks, ctr, face, (m, dt))
end

plot(ks, ctr)

ext: jl


text:
using WebSockets

const IP = "0.0.0.0"
const PORT = 8000
const KEY = "12345678"

function do_auth_res(data)
    @show data
    data == "AUTH:SUCCESS"
end

function do_handle_res(data)
    @show data
end

function send_message(message, ws_client)
    writeguarded(ws_client, "AUTH:$KEY")
    data, success = readguarded(ws_client)
    if !do_auth_res(data |> String)
        println("authentication fail!")
        return
    end
    writeguarded(ws_client, message)
    data, success = readguarded(ws_client)
    do_handle_res(data |> String)
end

send_message(message) = x->send_message(message, x)

MSG = "CNOT 0 1"
ws = WebSockets.open(send_message(MSG), "ws://$IP:$PORT")


text:
module Coping

using Codex
using DataFrames

export julliet2scores

ans_mapping = Dict(
    "Nooit" => 1,
    "Bijna nooit" => 2,
    "Soms" => 3,
    "Vaak" => 4,
    "Altijd" => 5
)

ans2num(ans)::Int = return ans_mapping[ans]

score(responses...) = sum(ans2num.(responses))

"""
    julliet2scores(df::DataFrame)

In 2021-07, returns 289×6 DataFrame.

For Julliet, 2018 and 2019 are the same.
2020 is a bit different, but it shouldn't matter for this code, see `diff`
For more information, see the e-mail with the subject `Julliet`.

This file is not tested in-depth, because the procedure is the same as delta.
"""
function julliet2scores(df::DataFrame)
    cols = [
        :id,
        :completed_at,
        [:v1, :v5, :v8, :v11, :v14] => ByRow(score) => :problem_focused,
        [:v3, :v7, :v10, :v13, :v4] => ByRow(score) => :emotion_focused,
        [:v2, :v6, :v9, :v15, :v12] => ByRow(score) => :seeking_support
    ]
    select!(df, cols...)
    disallowmissing!(df)
end

end # module


text:
using ProgressMeter, Plots, StatsBase, Statistics, LaTeXStrings, JLD

function TRW_CD(L,P,X0)
    surf = zeros(L+2); surf[X0] = 1.0
    TimePassed = 0
    while surf[1] + surf[end] < 0.99
        Altsurf = copy(surf)
        for i ∈ 2:L+1
            Altsurf[i+1] += surf[i]*(1-P)
            Altsurf[i-1] += surf[i]*P
            Altsurf[i] -= surf[i]
        end
        surf = Altsurf
        TimePassed += 1
    end
    return TimePassed
end

L = 40
PList = hcat(0:0.01:1)
X0List = hcat(1:L+2)

Data = zeros(length(X0List),length(PList))
progress = Progress(length(X0List)*length(PList)*runnum)

for i in 1:length(X0List)
    X0 = X0List[i]
    for j in 1:length(PList)
        P = PList[j]
        Data[i,j] = mean(TRW_CD(L,P,X0))
    end
end

save("../../Data/Q6/Q6-LT.jld", "Data", Data)
Data = load("../../Data/Q6/Q6-LT.jld")["Data"]

P1 = heatmap(0:0.01:1, 1:L+2,Data, xaxis=nothing,legend = nothing,title = L"Life\ time")
P2 = contour(0:0.01:1, 1:L+2,Data, fill = true ,legend = nothing,xlabel = L"P")
P3 = plot(P1,P2,layout = (2,1),ylabel = L"X_{init}",link = :y)
CBar = heatmap((0:1650).*ones(1651,1),xlabel ="",legend=:none, xticks=:none, yticks=(1:165:1651, string.(0:165:1651)))

plot([P3]...,CBar,layout = @layout[grid(1,1) a{0.04w}], size = (600,500))
savefig("../../Figs/Q6/Q6-LT.pdf")

X0List[21]

X0Plot = []
for X0 in [11,21,31]
    PLT = begin
        plot(PList,Data[X0,:,1], c = :black,linewidth = 3)
        plot!(PList,Data[X0,:,1], linewidth = 2, c = :purple, fill = (0, 0.5))
        plot!(title = L"X_{init} = %$(X0List[X0])",
            xlabel = L"P", ylabel = L"Avg\ Life\ Time", legend = :none)
    end
    push!(X0Plot, PLT)
end

X0P = plot(X0Plot...,layout = (1,3))

PPlot = []
for P in [26,51,76]
    PLT = begin
        plot(X0List,Data[:,P], c = :black, linewidth = 3)
        plot!(X0List,Data[:,P], linewidth = 2, c = :steelblue,  fill = (0, 0.5))
        plot!(title = L"P = %$(PList[P])",
            xlabel = L"X0", ylabel = L"Avg\ Life\ Time", legend = :none)
    end
    push!(PPlot, PLT)
end
PPlot[1]
PP = plot(PPlot...,layout = (1,3))

Pmain2 = plot(PP,X0P,layout = (2,1),size = (850,690), plot_title = L"Life\ time\ of\ a\ Random\ Walker\ (%$runnum \ runs)")
savefig(Pmain2,"../../Figs/Q6/Q6-scat.pdf")


text:
using ITensors,
      ITensorsGPU,
      Test
@testset "cuMPS Basics" begin

  N = 10
  sites = [Index(2,"Site") for n=1:N]
  psi = cuMPS(sites)
  @test length(psi) == N
  @test length(cuMPS()) == 0

  str = split(sprint(show, psi), '\n')
  @test str[1] == "MPS"
  @test length(str) == length(psi) + 2

  @test siteind(psi,2) == sites[2]
  @test hasind(psi[3],linkind(psi,2))
  @test hasind(psi[3],linkind(psi,3))

  psi[1] = cuITensor(sites[1])
  @test hasind(psi[1],sites[1])

  L = randomMPS(sites)
  K = cuMPS(L)
  @test all(ITensors.data(cpu(K)) .== ITensors.data(cpu(L)))

  @testset "cuproductMPS" begin
    @testset "vector of string input" begin
      sites = siteinds("S=1/2",N)
      state = fill("",N)
      for j=1:N
        state[j] = isodd(j) ? "Up" : "Dn"
      end
      psi  = productCuMPS(sites,state)
      for j=1:N
        sign = isodd(j) ? +1.0 : -1.0
        ops = cuITensor(op(sites,"Sz",j))
        psip = prime(psi[j], "Site")
        res = psi[j]*ops*dag(psip)
        @test res[] ≈ sign/2
      end
      @test_throws DimensionMismatch productCuMPS(sites, fill("", N - 1))
    end

    @testset "vector of int input" begin
      sites = siteinds("S=1/2",N)
      state = fill(0,N)
      for j=1:N
        state[j] = isodd(j) ? 1 : 2
      end
      psi = productCuMPS(sites,state)
      for j=1:N
        sign = isodd(j) ? +1.0 : -1.0
        ops = cuITensor(op(sites,"Sz",j))
        psip = prime(psi[j], "Site")
        @test (psi[j]*ops*dag(psip))[] ≈ sign/2
      end
    end

  end

  @testset "randomMPS" begin
    phi = randomCuMPS(sites)
    @test hasind(phi[1],sites[1])
    @test norm(phi[1])≈1.0
    @test hasind(phi[4],sites[4])
    @test norm(phi[4])≈1.0
  end

  @testset "inner different MPS" begin
    phi = randomMPS(sites)
    psi = randomMPS(sites)
    phipsi = dag(phi[1])*psi[1]
    for j = 2:N
      phipsi *= dag(phi[j])*psi[j]
    end
    @test phipsi[] ≈ inner(phi,psi)
    phi = randomCuMPS(sites)
    psi = randomCuMPS(sites)
    cphi = MPS([cpu(phi[i]) for i in 1:length(phi)])
    cpsi = MPS([cpu(psi[i]) for i in 1:length(psi)])
    phipsi = dag(phi[1])*psi[1]
    cphipsi = dag(cphi[1])*cpsi[1]
    for j = 2:N
      phipsi  *= dag(phi[j])*psi[j]
      cphipsi *= dag(cphi[j])*cpsi[j]
    end
    @test cpu(phipsi)[] ≈ cphipsi[]
    @test cpu(phipsi)[] ≈ inner(cphi,cpsi)
    @test cpu(phipsi)[] ≈ inner(phi,psi)
    phipsi = dag(phi[1])*psi[1]
    for j = 2:N
      phipsi = phipsi*dag(phi[j])*psi[j]
    end
    @test cpu(phipsi)[] ≈ inner(phi,psi)

    badsites = [Index(2) for n=1:N+1]
    badpsi = randomCuMPS(badsites)
    @test_throws DimensionMismatch inner(phi,badpsi)
  end

  @testset "inner same MPS" begin
    psi    = randomMPS(sites)
    psidag = dag(deepcopy(psi))
    ITensors.prime_linkinds!(psidag)
    psipsi = psidag[1]*psi[1]
    for j = 2:N
      psipsi = psipsi * psidag[j]*psi[j]
    end
    @test psipsi[] ≈ inner(psi,psi)
    psi    = randomCuMPS(sites)
    psidag = dag(deepcopy(psi))
    ITensors.prime_linkinds!(psidag)
    psipsi = psidag[1]*psi[1]
    for j = 2:N
      psipsi = psipsi *psidag[j]*psi[j]
    end
    @test psipsi[] ≈ inner(psi,psi)
  end

  @testset "add MPS" begin
    psi      = randomMPS(sites)
    phi      = deepcopy(psi)
    xi       = add(psi, phi)
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
    psi      = randomCuMPS(sites)
    phi      = deepcopy(psi)
    xi       = add(psi, phi)
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
  end

  sites = siteinds("S=1/2",N)
  psi = cuMPS(sites)
  @test length(psi) == N # just make sure this works
  @test length(siteinds(psi)) == N

  psi = randomCuMPS(sites)
  orthogonalize!(psi, N-1)
  @test ITensors.leftlim(psi) == N-2
  @test ITensors.rightlim(psi) == N
  orthogonalize!(psi, 2)
  @test ITensors.leftlim(psi) == 1
  @test ITensors.rightlim(psi) == 3
  psi = randomCuMPS(sites)
  psi.rlim = N+1 # do this to test qr from rightmost tensor
  orthogonalize!(psi, div(N, 2))
  @test ITensors.leftlim(psi) == div(N, 2) - 1
  @test ITensors.rightlim(psi) == div(N, 2) + 1

  #@test_throws ErrorException linkind(MPS(N, fill(cuITensor(), N), 0, N + 1), 1)

  @testset "replacebond!" begin
  # make sure factorization preserves the bond index tags
    psi = randomCuMPS(sites)
    phi = psi[1]*psi[2]
    bondindtags = tags(linkind(psi,1))
    replacebond!(psi,1,phi)
    @test tags(linkind(psi,1)) == bondindtags

    # check that replaceBond! updates llim_ and rlim_ properly
    orthogonalize!(psi,5)
    phi = psi[5]*psi[6]
    replacebond!(psi,5,phi, ortho="left")
    @test ITensors.leftlim(psi)==5
    @test ITensors.rightlim(psi)==7

    phi = psi[5]*psi[6]
    replacebond!(psi,5,phi,ortho="right")
    @test ITensors.leftlim(psi)==4
    @test ITensors.rightlim(psi)==6

    psi.llim = 3
    psi.rlim = 7
    phi = psi[5]*psi[6]
    replacebond!(psi,5,phi,ortho="left")
    @test ITensors.leftlim(psi)==3
    @test ITensors.rightlim(psi)==7
  end

end

# Helper function for making MPS
function basicRandomCuMPS(N::Int;dim=4)
  sites  = [Index(2,"Site") for n=1:N]
  M      = MPS(sites)
  links  = [Index(dim,"n=$(n-1),Link") for n=1:N+1]
  M[1]   = randomCuITensor(sites[1],links[2])
  for n=2:N-1
    M[n] = randomCuITensor(links[n],sites[n],links[n+1])
  end
  M[N]   = randomCuITensor(links[N],sites[N])
  M[1]  /= sqrt(inner(M,M))
  return M
end
@testset "MPS gauging and truncation" begin

  N = 30

  @testset "orthogonalize! method" begin
    c = 12
    M = basicRandomCuMPS(N)
    orthogonalize!(M,c)

    @test ITensors.leftlim(M) == c-1
    @test ITensors.rightlim(M) == c+1

    # Test for left-orthogonality
    L = M[1]*prime(M[1],"Link")
    l = linkind(M,1)
    @test cpu(L) ≈ delta(l,l') rtol=1E-12
    for j=2:c-1
      L = L*M[j]*prime(M[j],"Link")
      l = linkind(M,j)
      @test cpu(L) ≈ delta(l,l') rtol=1E-12
    end

    # Test for right-orthogonality
    R = M[N]*prime(M[N],"Link")
    r = linkind(M,N-1)
    @test cpu(R) ≈ delta(r,r') rtol=1E-12
    for j in reverse(c+1:N-1)
      R = R*M[j]*prime(M[j],"Link")
      r = linkind(M,j-1)
      @test cpu(R) ≈ delta(r,r') rtol=1E-12
    end

    @test norm(M[c]) ≈ 1.0
  end

  @testset "truncate! method" begin
    M  = basicRandomCuMPS(N;dim=10)
    M0 = copy(M)
    truncate!(M;maxdim=5)
    @test ITensors.rightlim(M) == 2
    # Test for right-orthogonality
    R = M[N]*prime(M[N],"Link")
    r = linkind(M,N-1)
    @test cpu(R) ≈ delta(r,r') rtol=1E-12
    for j in reverse(2:N-1)
      R = R*M[j]*prime(M[j],"Link")
      r = linkind(M,j-1)
      @test cpu(R) ≈ delta(r,r') rtol=1E-12
    end
    @test inner(M0, M) > 0.1
  end

end

#=@testset "Other MPS methods" begin

  @testset "sample! method" begin
    N = 10
    sites = [Index(3,"Site,n=$n") for n=1:N]
    psi = makeRandomCuMPS(sites,chi=3)
    nrm2 = inner(psi,psi)
    psi[1] *= (1.0/sqrt(nrm2))

    s = sample!(psi)

    @test length(s) == N
    for n=1:N
      @test 1 <= s[n] <= 3
    end

    # Throws becase not orthogonalized to site 1:
    orthogonalize!(psi,3)
    @test_throws ErrorException sample(psi)

    # Throws becase not normalized
    orthogonalize!(psi,1)
    psi[1] *= (5.0/norm(psi[1]))
    @test_throws ErrorException sample(psi)

    # Works when ortho & normalized:
    orthogonalize!(psi,1)
    psi[1] *= (1.0/norm(psi[1]))
    s = sample(psi)
    @test length(s) == N
  end

end=#


text:
export Nearest, Previous, Next

abstract type ConstantInterpType end
struct Nearest <: ConstantInterpType end
struct Previous <: ConstantInterpType end
struct Next <: ConstantInterpType end

struct Constant{T<:ConstantInterpType} <: Degree{0}
    Constant() = new{Nearest}()
    Constant{Nearest}() = new{Nearest}()
    Constant{Previous}() = new{Previous}()
    Constant{Next}() = new{Next}()
end

"""
Constant b-splines are *nearest-neighbor* interpolations, and effectively
return `A[round(Int,x)]` when interpolating.
"""
Constant

function positions(c::Constant{Previous}, ax, x)  # discontinuity occurs at integer locations
    xm = floorbounds(x, ax)
    δx = x - xm
    fast_trunc(Int, xm), δx
end
function positions(c::Constant{Next}, ax, x)  # discontinuity occurs at integer locations
    xm = ceilbounds(x, ax)
    δx = x - xm
    fast_trunc(Int, xm), δx
end
function positions(c::Constant{Nearest}, ax, x)  # discontinuity occurs at half-integer locations
    xm = roundbounds(x, ax)
    δx = x - xm
    fast_trunc(Int, xm), δx
end

value_weights(::Constant, δx) = (1,)
gradient_weights(::Constant, δx) = (0,)
hessian_weights(::Constant, δx) = (0,)

padded_axis(ax::AbstractUnitRange, ::BSpline{<:Constant}) = ax


text:
using LinearAlgebra
using JuMP

function test(n::Int)
    model = Model()
    a = rand(n, n)
    @variable(model, x[1:n, 1:n])
    b = rand(n, n, n)
    @variable(model, y[1:n, 1:n, 1:n])
    c = rand(n)
    @variable(model, z[1:n])

    #initialize
    @constraint(model, sum(c[i] * z[i] for i in 1:n) <= 0)

    #Vector
    elapsed = @elapsed @constraint(model, sum(c[i] * z[i] for i in 1:n) <= 1)
    println("Vector with sum(): $(elapsed)")

    elapsed = @elapsed @constraint(model, dot(c, z) <= 1)
    println("Vector with vecdot() : $(elapsed)")

    #2D Matrix
    elapsed = @elapsed @constraint(
        model,
        sum(a[i, j] * x[i, j] for i in 1:n, j in 1:n) <= 1
    )
    println("2D Matrix with sum(): $(elapsed)")

    elapsed = @elapsed @constraint(model, dot(a, x) <= 1)
    println("2D Matrix with bigvecdot(): $(elapsed)")

    #3D Matrix
    elapsed = @elapsed @constraint(
        model,
        sum(b[i, j, k] * y[i, j, k] for i in 1:n, j in 1:n, k in 1:n) <= 1
    )
    println("3D Matrix with sum(): $(elapsed)")

    elapsed = @elapsed @constraint(model, dot(b, y) <= 1)
    println("3D Matrix with vecdot(): $(elapsed)")
    return 0
end

for n in [10, 50, 100, 200, 300]
    println("n = $n")
    test(n)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function get_test_trajdata(roadway::Roadway)
    trajdata = Trajdata(0.1)

    trajdata.defs[1] = VehicleDef(AgentClass.CAR, 5.0, 3.0)
    trajdata.defs[2] = VehicleDef(AgentClass.CAR, 5.0, 3.0)

    push!(trajdata.states, RecordState{VehicleState, Int64}(VehicleState(VecSE2(0.0,0.0,0.0), roadway, 10.0), 1)) # car 1, frame 1
    push!(trajdata.states, RecordState{VehicleState, Int64}(VehicleState(VecSE2(3.0,0.0,0.0), roadway, 20.0), 2)) # car 2, frame 1
    push!(trajdata.states, RecordState{VehicleState, Int64}(VehicleState(VecSE2(1.0,0.0,0.0), roadway, 10.0), 1)) # car 1, frame 2
    push!(trajdata.states, RecordState{VehicleState, Int64}(VehicleState(VecSE2(5.0,0.0,0.0), roadway, 20.0), 2)) # car 2, frame 2

    push!(trajdata.frames, RecordFrame(1,2))
    push!(trajdata.frames, RecordFrame(3,4))

    trajdata
end

@testset "1d state" begin
    s = State1D(0.0, 0.0)
    path, io = mktemp()
    write(io, MIME"text/plain"(), s)
    close(io)
    io = open(path)
    s2 = read(io, MIME"text/plain"(), State1D)
    close(io)
    @test s == s2

    veh = Vehicle1D(s, VehicleDef(), 1)
    scene = Scene1D()
    push!(scene, veh)
    scene2 = Scene1D([veh])
    @test scene[1].state.s == 0.0
    @test first(scene2.entities) == first(scene.entities)
    @test scene2.n == scene.n == 1
    @test get_center(veh) == 0.0
    @test get_footpoint(veh) == 0.0
    @test get_front(veh) == veh.def.length/2
    @test get_rear(veh) == - veh.def.length/2

end

@testset "VehicleState" begin
    s = VehicleState(VecSE2(0.0,0.0,0.0), Frenet(NULL_ROADINDEX, 0.0, 0.0, 0.0), 10.0)
    @test isapprox(velf(s).s, 10.0)
    @test isapprox(velf(s).t,  0.0)
    show(IOBuffer(), s.posF)
    show(IOBuffer(), s)

    s = VehicleState(VecSE2(0.0,0.0,0.0), Frenet(NULL_ROADINDEX, 0.0, 0.0, 0.1), 10.0)
    @test isapprox(get_vel_s(s), 10.0*cos(0.1))
    @test isapprox(get_vel_t(s), 10.0*sin(0.1))
    @test isapprox(velf(s).s, 10.0*cos(0.1))
    @test isapprox(velf(s).t, 10.0*sin(0.1))
    @test isapprox(velg(s).x, 10.0)
    @test isapprox(velg(s).y, 0.0)
    @test isapprox(vel(s), 10.0)

    vehdef = VehicleDef(AgentClass.CAR, 5.0, 3.0)
    veh = Vehicle(s, vehdef, 1)
    @test isapprox(get_footpoint(veh), VecSE2(0.0,0.0,0.0))
    show(IOBuffer(), vehdef)
    show(IOBuffer(), veh)

    ri = RoadIndex(CurveIndex(1,0.1), LaneTag(1,2))
    @test isapprox(Frenet(ri, 0.0, 0.0, 0.1), Frenet(ri, 0.0, 0.0, 0.1))
    @test VehicleState(VecSE2(0.1,0.2,0.3), 1.0) == VehicleState(VecSE2(0.1,0.2,0.3), NULL_FRENET, 1.0)

    @test get_front(veh).x == vehdef.length/2
    @test get_rear(veh).x == -vehdef.length/2

    roadway = gen_straight_roadway(3, 100.0)
    veh = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 0.0)
    veh = Vehicle(veh, VehicleDef(), 1)
    @test get_lane(roadway, veh).tag == LaneTag(1,1)
    @test get_lane(roadway, veh).tag == get_lane(roadway, veh.state).tag
    veh = VehicleState(VecSE2(0.0, 3.0, 0.0), roadway, 0.0)
    @test get_lane(roadway, veh).tag == LaneTag(1,2)
    veh = VehicleState(VecSE2(0.0, 6.0, 0.0), roadway, 0.0)
    @test get_lane(roadway, veh).tag == LaneTag(1,3)
end

@testset "scene" begin
    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)

    s1 = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 0.0)
    s2 = VehicleState(VecSE2(5.0, 0.0, 0.0), roadway, 0.0)
    s3 = lerp(s1, s2, 0.5, roadway)
    @test s3.posG.x == 2.5
    @test s3.posF.s == 2.5

    vehstate = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 0.0)
    vehstate1 = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway[LaneTag(1,1)], roadway, 0.0)
    @test vehstate1 == vehstate
    veh = Vehicle(vehstate, VehicleDef(), 1)
    scene1 = Scene()
    push!(scene1, veh)
    scene2 = Scene([veh])
    @test first(scene1.entities) == first(scene2.entities)
    @test scene1.n == scene2.n

    io = IOBuffer()
    show(io, scene1)
    close(io)

    veh2 = Entity(vehstate, BicycleModel(VehicleDef()), 1)
    veh3 = convert(Vehicle, veh2)
    @test veh3 == veh

    rec = SceneRecord(1, 0.5)
    rec = SceneRecord(1, 0.5, 10)

    io = IOBuffer()
    show(io, rec)
    close(io)
    scene = Scene()
    get!(scene, trajdata, 1)
    @test length(scene) == 2
    for (i,veh) in enumerate(scene)
        @test scene[i].state == get_state(trajdata, i, 1)
        @test scene[i].def == get_def(trajdata, i)
    end

    scene2 = Scene(deepcopy(scene.entities), 2)
    @test length(scene2) == 2
    for (i,veh) in enumerate(scene2)
        @test scene2[i].state == get_state(trajdata, i, 1)
        @test scene2[i].def == get_def(trajdata, i)
    end

    @test get_by_id(scene, 1) == scene[1]

    empty!(scene2)
    @test length(scene2) == 0

    copyto!(scene2, scene)
    @test length(scene2) == 2
    for (i,veh) in enumerate(scene2)
        @test scene2[i].state == get_state(trajdata, i, 1)
        @test scene2[i].def == get_def(trajdata, i)
    end

    delete!(scene2, scene2[1])
    @test length(scene2) == 1
    @test scene2[1].state == get_state(trajdata, 2, 1)
    @test scene2[1].def == get_def(trajdata, 2)
    scene2[1] = deepcopy(scene[1])
    @test scene2[1].state == get_state(trajdata, 1, 1)
    @test scene2[1].def == get_def(trajdata, 1)

    @test findfirst(1, scene) == 1
    @test findfirst(2, scene) == 2
    @test get_first_available_id(scene) == 3

    @test in(1, scene)
    @test in(2, scene)
    @test !in(3, scene)

    veh = scene[2]
    @test veh.state == get_state(trajdata, 2, 1)
    @test veh.def == get_def(trajdata, 2)

    push!(scene, get_state(trajdata, 1, 1))
end

@testset "trajdata" begin
    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)

    @test nframes(trajdata) == 2
    @test !frame_inbounds(trajdata, 0)
    @test frame_inbounds(trajdata, 1)
    @test frame_inbounds(trajdata, 2)
    @test !frame_inbounds(trajdata, 3)

    @test n_objects_in_frame(trajdata, 1) == 2
    @test n_objects_in_frame(trajdata, 2) == 2

    @test nth_id(trajdata, 1) == 1
    @test nth_id(trajdata, 1, 2) == 2
    @test nth_id(trajdata, 2, 1) == 1
    @test nth_id(trajdata, 2, 2) == 2

    @test findfirst_frame_with_id(trajdata, 1) == 1
    @test findfirst_frame_with_id(trajdata, 2) == 1
    @test findfirst_frame_with_id(trajdata, -1) == nothing
    @test findlast_frame_with_id(trajdata, 1) == 2
    @test findlast_frame_with_id(trajdata, 2) == 2
    @test findlast_frame_with_id(trajdata, -1) == nothing

    @test sort!(get_ids(trajdata)) == [1,2]

    @test in(1, trajdata, 1)
    @test in(1, trajdata, 2)
    @test in(2, trajdata, 1)
    @test in(2, trajdata, 2)
    @test !in(3, trajdata, 1)

    @test isapprox(get_time(trajdata, 1), 0.0)
    @test isapprox(get_time(trajdata, 2), 0.1)

    @test isapprox(get_elapsed_time(trajdata, 1, 2),  0.1)
    @test isapprox(get_elapsed_time(trajdata, 2, 1), -0.1)

    @test get_timestep(trajdata) == 0.1

    veh = get(trajdata, 1, 1)
    @test veh.state == VehicleState(VecSE2(0.0,0.0,0.0), roadway, 10.0)
    @test_throws ArgumentError get(trajdata, 10, 1)
    @test_throws BoundsError get(trajdata, 1, 10)

    let
        iter = ListRecordStateByIdIterator(trajdata, 1)
        items = collect(iter) # list of (frame_index, state)
        @test length(items) == 2
        @test items[1][1] == 1
        @test items[1][2] == get_state(trajdata, 1, 1)
        @test items[2][1] == 2
        @test items[2][2] == get_state(trajdata, 1, 2)

        iter = ListRecordStateByIdIterator(trajdata, 2)
        items = collect(iter)
        @test length(items) == 2
        @test items[1][1] == 1
        @test items[1][2] == get_state(trajdata, 2, 1)
        @test items[2][1] == 2
        @test items[2][2] == get_state(trajdata, 2, 2)
    end

    path, io = mktemp()
    write(io, MIME"text/plain"(), trajdata)
    close(io)

    io = open(path)
    trajdata2 = read(io, MIME"text/plain"(), Trajdata)
    close(io)
    rm(path)

    @test nframes(trajdata2) == nframes(trajdata)
    for i in 1 : nframes(trajdata2)
        @test n_objects_in_frame(trajdata2, i) == n_objects_in_frame(trajdata, i)
        for j in 1 : n_objects_in_frame(trajdata, i)
            veh1 = get(trajdata, j, i)
            veh2 = get(trajdata2, j, i)
            @test veh1.id == veh2.id
            @test veh1.def.class == veh2.def.class
            @test isapprox(veh1.def.length, veh2.def.length)
            @test isapprox(veh1.def.width, veh2.def.width)

            @test isapprox(veh1.state.v, veh2.state.v)
            @test isapprox(veh1.state.posG, veh2.state.posG, atol=1e-3)
            @test isapprox(veh1.state.posF.s, veh2.state.posF.s, atol=1e-3)
            @test isapprox(veh1.state.posF.t, veh2.state.posF.t, atol=1e-3)
            @test isapprox(veh1.state.posF.ϕ, veh2.state.posF.ϕ, atol=1e-6)
            @test veh1.state.posF.roadind.tag == veh2.state.posF.roadind.tag
            @test veh1.state.posF.roadind.ind.i == veh2.state.posF.roadind.ind.i
            @test isapprox(veh1.state.posF.roadind.ind.t, veh2.state.posF.roadind.ind.t, atol=1e-3)
        end
    end

    trajdata3 = get_subinterval(trajdata2, 1, nframes(trajdata2))
    @test nframes(trajdata3) == nframes(trajdata2)
    for i in 1 : nframes(trajdata3)
        @test n_objects_in_frame(trajdata3, i) == n_objects_in_frame(trajdata2, i)
        for j in 1 : n_objects_in_frame(trajdata2, i)
            veh1 = get(trajdata2, j, i)
            veh2 = get(trajdata3, j, i)
            @test veh1.id == veh2.id
            @test veh1.def.class == veh2.def.class
            @test isapprox(veh1.def.length, veh2.def.length)
            @test isapprox(veh1.def.width, veh2.def.width)

            @test isapprox(veh1.state.v, veh2.state.v)
            @test isapprox(veh1.state.posG, veh2.state.posG, atol=1e-3)
            @test isapprox(veh1.state.posF.s, veh2.state.posF.s, atol=1e-3)
            @test isapprox(veh1.state.posF.t, veh2.state.posF.t, atol=1e-3)
            @test isapprox(veh1.state.posF.ϕ, veh2.state.posF.ϕ, atol=1e-6)
            @test veh1.state.posF.roadind.tag == veh2.state.posF.roadind.tag
            @test veh1.state.posF.roadind.ind.i == veh2.state.posF.roadind.ind.i
            @test isapprox(veh1.state.posF.roadind.ind.t, veh2.state.posF.roadind.ind.t, atol=1e-3)
        end
    end

    trajdata3 = get_subinterval(trajdata2, 1, 1)
    @test nframes(trajdata3) == 1
    let
        i = 1
        @test n_objects_in_frame(trajdata3, i) == n_objects_in_frame(trajdata2, i)
        for j in 1 : n_objects_in_frame(trajdata2, i)
            veh1 = get(trajdata2, j, i)
            veh2 = get(trajdata3, j, i)
            @test veh1.id == veh2.id
            @test veh1.def.class == veh2.def.class
            @test isapprox(veh1.def.length, veh2.def.length)
            @test isapprox(veh1.def.width, veh2.def.width)

            @test isapprox(veh1.state.v, veh2.state.v)
            @test isapprox(veh1.state.posG, veh2.state.posG, atol=1e-3)
            @test isapprox(veh1.state.posF.s, veh2.state.posF.s, atol=1e-3)
            @test isapprox(veh1.state.posF.t, veh2.state.posF.t, atol=1e-3)
            @test isapprox(veh1.state.posF.ϕ, veh2.state.posF.ϕ, atol=1e-6)
            @test veh1.state.posF.roadind.tag == veh2.state.posF.roadind.tag
            @test veh1.state.posF.roadind.ind.i == veh2.state.posF.roadind.ind.i
            @test isapprox(veh1.state.posF.roadind.ind.t, veh2.state.posF.roadind.ind.t, atol=1e-3)
        end
    end
end

@testset "SceneRecord" begin
    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)

    Δt = 0.1
    rec = SceneRecord(5, Δt)
    @test capacity(rec) == 5
    @test nframes(rec) == 0
    @test !pastframe_inbounds(rec, 0)
    @test !pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)

    scene = get!(Scene(), trajdata, 1)
    update!(rec, scene)
    @test nframes(rec) == 1
    @test pastframe_inbounds(rec, 0)
    @test !pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)
    @test isapprox(get_elapsed_time(rec, 0), 0)
    @test rec[0][1].state == get_state(trajdata, 1, 1)
    @test rec[0][1].def == get_def(trajdata, 1)
    @test rec[0][2].state == get_state(trajdata, 2, 1)
    @test rec[0][2].def == get_def(trajdata, 2)
    show(IOBuffer(), rec)

    get!(scene, trajdata, 2)
    update!(rec, scene)
    @test nframes(rec) == 2
    @test pastframe_inbounds(rec, 0)
    @test pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)
    @test isapprox(get_elapsed_time(rec,  0),  0)
    @test isapprox(get_elapsed_time(rec, -1), Δt)
    @test isapprox(get_elapsed_time(rec, -1, 0), Δt)
    @test rec[0][1].state == get_state(trajdata, 1, 2)
    @test rec[0][1].def == get_def(trajdata, 1)
    @test rec[0][2].state == get_state(trajdata, 2, 2)
    @test rec[0][2].def == get_def(trajdata, 2)
    @test rec[-1][1].state == get_state(trajdata, 1, 1)
    @test rec[-1][1].def == get_def(trajdata, 1)
    @test rec[-1][2].state == get_state(trajdata, 2, 1)
    @test rec[-1][2].def == get_def(trajdata, 2)

    scene2 = get!(Scene(), rec)
    @test scene2[1].state == get_state(trajdata, 1, 2)
    @test scene2[1].def == get_def(trajdata, 1)
    @test scene2[2].state == get_state(trajdata, 2, 2)
    @test scene2[2].def == get_def(trajdata, 2)

    get!(scene2, rec, -1)
    @test scene2[1].state == get_state(trajdata, 1, 1)
    @test scene2[1].def == get_def(trajdata, 1)
    @test scene2[2].state == get_state(trajdata, 2, 1)
    @test scene2[2].def == get_def(trajdata, 2)

    empty!(rec)
    @test nframes(rec) == 0

    test_veh_state = VehicleState(VecSE2(7.0,7.0,2.0), roadway, 10.0)
    test_veh_def = VehicleDef(AgentClass.CAR, 5.0, 3.0)
    test_veh = Vehicle(test_veh_state, test_veh_def, 999)
    rec[-1][1] = test_veh
    @test rec[-1][1].state == test_veh_state
end


text:
export d2u, j2md, md2j, parsetimewin, timestamp, u2d, timespec

function tstr(t::DateTime)
  Y, M, D, h, m, s, μ = year(t), month(t), day(t), hour(t), minute(t), second(t), millisecond(t)
  Y = lpad(Y, 4, "0")
  M = lpad(M, 2, "0")
  D = lpad(D, 2, "0")
  h = lpad(h, 2, "0")
  m = lpad(m, 2, "0")
  s = lpad(s, 2, "0")
  μ = lpad(μ, 3, "0")
  return string(Y, "-", M, "-", D, "T", h, ":", m, ":", s, ".", μ)
end

"Alias to Dates.unix2datetime"
u2d(k::Real) = Dates.unix2datetime(k)
"Alias to Dates.datetime2unix"
d2u(k::DateTime) = Dates.datetime2unix(k)

@doc """
    timestamp()

Return current time formatted YYYY-mm-ddTHH:MM:SS.nnn.
""" timestamp
timestamp() = tstr(Dates.unix2datetime(time()))
timestamp(t::DateTime) = tstr(t)
timestamp(t::Real) = tstr(u2d(t))
timestamp(t::String) = tstr(Dates.DateTime(t))
tnote(s::String) = string(timestamp(), ": ", s)

"""
  m,d = j2md(y,j)

Convert Julian day j of year y to month m, day d
"""
function j2md(y::T, j::T) where T<:Integer
  if T != Int32
    y = Int32(y)
    j = Int32(j)
  end
  z = zero(Int32)
  o = one(Int32)
  m = z
  d = o
  if j > Int32(31)
    leapyear = ((j > 59) && ((y % Int32(400) == z) || (y % Int32(4) == z && y % Int32(100) != z)))
    while j > z
      d = j
      m += o
      j -= days_per_month[m]
      if leapyear && m == 2
        j -= o
      end
    end
  else
    m = o
    d = j
  end
  return m,d
end

"""
  j = md2j(y,m,d)

Convert month `m`, day `d` of year `y` to Julian day (day of year)
"""
function md2j(y::T, m::T, d::T) where T<:Integer
  j = zero(Int32)
  if T != Int32
    y = Int32(y)
    m = Int32(m)
    d = Int32(d)
  end
  z = zero(Int32)
  i = one(Int8)
  while i < m
    j = j+getindex(days_per_month, i)
    i = i+1
  end
  j = j+d
  if m > 2 && ((y % Int32(400) == z) ||
               (y % Int32(4)   == z &&
                y % Int32(100) != z))
    j = j+1
  end
  return T(j)
end
md2j(y::AbstractString, m::AbstractString, d::AbstractString) = md2j(parse(Int, y), parse(Int, m), parse(Int, d))

@doc """
# Time Specification

Most functions that allow time specifiation use two reserved key words to track
time:
* `s`: Start (begin) time
* `t`: Termination (end) time

The values passed to keywords `s` and `t` can be real numbers, DateTime objects,
or ASCII strings. Strings must follow the format "yyyy-mm-ddTHH:MM:SS.nnn", e.g.
`s="2016-03-23T11:17:00.333"`. Exact behavior depends on the types of s and t:

| **s** | **t** | **Behavior**                         |
|:------|:------|:-------------------------------------|
| DT    | DT    | Sort                                 |
| R     | DT    | Add `s` seconds to `t`, then sort    |
| DT    | R     | Add `t` seconds to `s`, then sort    |
| S     | R, DT | Convert `s` → DateTime, then sort    |
| R, DT | S     | Convert `t` → DateTime, then sort    |
| R     | R     | `s` and `t` are seconds from `now()` |

(above, R = Real, DT = DateTime, S = String, I = Integer)

## Relative Timekeeping
Numeric time values are interpreted *relative to the start of the current minute*.
Thus, if `-s` or `-t` is 0, the data request begins (or ends) at the start of
the minute in which the request is submitted.
""" timespec
timespec() = nothing

"""
d0, d1 = parsetimewin(s, t)

Convert times `s` and `t` to strings and sorts s.t. d0 < d1.

See also: ?timespec
"""
function parsetimewin(s::DateTime, t::DateTime)
  if s < t
    return (string(s), string(t))
  else
    return (string(t), string(s))
  end
end
parsetimewin(s::DateTime, t::String) = parsetimewin(s, DateTime(t))
parsetimewin(s::DateTime, t::Real) = parsetimewin(s, u2d(d2u(s)+t))
parsetimewin(s::Real, t::DateTime) = parsetimewin(t, u2d(d2u(t)+s))
parsetimewin(s::String, t::Union{Real,DateTime}) = parsetimewin(DateTime(s), t)
parsetimewin(s::Union{Real,DateTime}, t::String) = parsetimewin(s, DateTime(t))
parsetimewin(s::String, t::String) = parsetimewin(DateTime(s), DateTime(t))
parsetimewin(s::Real, t::Real) = parsetimewin(u2d(60*floor(Int, time()/60) + s), u2d(60*floor(Int, time()/60) + t))

# =========================================================
# Not for export
function y2μs(y::T) where T<:Integer
  y = Int64(y)-1
  return 86400000000 * (y*365 + div(y,4) - div(y,100) + div(y,400)) - 62135596800000000
end

# ts = round(Int64, d2u(DateTime(iv[1], m, d, iv[3], iv[4], iv[5], iv[6]))*sμ
mktime(y::T, j::T, h::T, m::T, s::T, μ::T) where T<:Integer = (y2μs(y) +
          Int64(j-one(T))*86400000000 +
          Int64(h)*3600000000 +
          Int64(m)*60000000 +
          Int64(s)*1000000 +
          Int64(μ))
mktime(t::Array{T,1}) where T<:Integer =(y2μs(t[1]) +
          Int64(t[2]-one(T))*86400000000 +
          Int64(t[3])*3600000000 +
          Int64(t[4])*60000000 +
          Int64(t[5])*1000000 +
          Int64(t[6]))

# convert a formatted time string to integer μs from the Unix epoch
function tstr2int(s::String)
  str = split(s, ".", limit=2)
  if length(str) < 2
    μ = 0
  else
    μ = parse(Int64, rpad(str[2], 6, '0'))
  end
  return DateTime(str[1]).instant.periods.value*1000 - dtconst + μ
end

# convert a time in integer μs (measured from the Unix epoch) to a string
function int2tstr(t::Int64)
  dt = unix2datetime(div(t, 1000000))
  v = 1000*getfield(getfield(getfield(dt, :instant), :periods), :value)
  r = string(t - v + dtconst)
  s = string(dt) * "." * lpad(r, 6, '0')
  return s
end

# =========================================================
# Time windowing functions
"""
    endtime(t::Array{Int64,2}, Δ::Int64)

Compute the time of the last sample in SeisIO time matrix `t` sampled at
(integer μs sampling interval) Δ. Result is measured from the Unix epoch in
integer μs.
"""
function endtime(t::Array{Int64,2}, Δ::Int64)
  if isempty(t)
    t_end = 0
  else
    L = size(t,1)
    t_end = (t[L,1]-1)*Δ
    if L > 2
      t_end += getindex(sum(t, dims=1),2)
    else
      t_end += t[1,2]
    end
    # t_end = getindex(sum(t, dims=1),2) + (t[L,1]-1)*Δ
  end
  return t_end
end
function endtime(t::Array{Int64,2}, fs::Float64)
  if fs == 0.0
    return t[size(t,1), 2]
  else
    return endtime(t, round(Int64, 1.0/(fs*μs)))
  end
end

function t_expand(t::Array{Int64,2}, fs::Float64)
  fs == 0.0 && return t[:,2]
  t[end,1] == 1 && return [t[1,2]]
  dt = round(Int64, 1.0/(fs*μs))
  tt = dt.*ones(Int64, t[end,1])
  tt[1] -= dt
  for i = 1:size(t,1)
    tt[t[i,1]] += t[i,2]
  end
  cumsum!(tt, tt)
  return tt
end

function t_collapse(tt::Array{Int64,1}, fs::Float64)
  if fs == 0.0
    t = hcat(collect(1:1:length(tt)), tt)
  else
    dt = round(Int64, 1.0/(fs*μs))
    ts = Array{Int64,1}([dt; diff(tt)::Array{Int64,1}])
    L = length(tt)
    i = findall(ts .!= dt)
    t = Array{Int64,2}([[1 tt[1]];[i ts[i].-dt]])
    if isempty(i) || i[end] != L
      t = vcat(t, hcat(L,0))
    end
  end
  return t
end

function x_inds(t::Array{Int64,2})
  nt = size(t, 1)-1
  inds = zeros(Int64, nt, 2)
  for i in 1:nt
    inds[i,1] = t[i,1]
    inds[i,2] = t[i+1,1] - (i == nt ? 0 : 1)
  end
  return inds
end

function t_win(T::Array{Int64,2}, Δ::Int64)
  isempty(T) && return(T)
  n = size(T,1)-1
  if T[n+1,2] != 0
    T = vcat(T, [T[n+1,1] 0])
    n += 1
  end
  w0 = -(Δ)
  W = Array{Int64,2}(undef,n,2)
  for i = 1:n
    W[i,1] = T[i,2] + w0 + Δ
    W[i,2] = W[i,1] + Δ*(T[i+1,1]-T[i,1]-1)
    w0 = W[i,2]
  end
  W[n,2] += Δ
  return W
end
t_win(T::Array{Int64,2}, fs::Float64) = t_win(T, round(Int64, 1000000.0/fs))

function w_time(W::Array{Int64,2}, Δ::Int64)
  n = size(W,1)+1
  T = Array{Int64,2}(undef,n,2)
  T[1,1] = Int64(1)
  T[1,2] = W[1,1]
  for i = 2:n-1
    T[i,1] = T[i-1,1] + div(W[i-1,2]-W[i-1,1], Δ) + 1
    T[i,2] = W[i,1] - W[i-1,2] - Δ
  end
  T[n,1] = T[n-1,1] + div(W[n-1,2]-W[n-1,1], Δ)
  T[n,2] = 0
  if T[n,1] == T[n-1,1]
    T = T[1:n-1,:]
  end
  return T
end
w_time(W::Array{Int64,2}, fs::Float64) = w_time(W, round(Int64, 1000000.0/fs))

function unpack_u8(v::UInt8)
  a = signed(div(v,0x10))*Int8(10)
  b = signed(rem(v,0x10))
  return Int64(a+b)
end

function datehex2μs!(a::Array{Int64,1}, datehex::Array{UInt8,1})
  a[1] = 100*unpack_u8(getindex(datehex, 1)) + unpack_u8(getindex(datehex, 2))
  a[2] = md2j(getindex(a,1), unpack_u8(getindex(datehex, 3)), unpack_u8(getindex(datehex, 4))) - 1
  setindex!(a, y2μs(getindex(a,1)), 1)
  a[3] = unpack_u8(getindex(datehex, 5))
  a[4] = unpack_u8(getindex(datehex, 6))
  a[5] = unpack_u8(getindex(datehex, 7))
  a[6] = unpack_u8(getindex(datehex, 8))
  return a[1] + a[2]*86400000000 + a[3]*3600000000 + a[4]*60000000 + a[5]*1000000 + a[6]*10000
end

function tx_float(t::Array{Int64,2}, fs::Float64)
  fs == 0.0 && return map(Float64, t[:,2])
  t[end,1] == 1 && return Float64[t[1,2]]
  Nt = t[end,1]
  dt = 1.0/fs
  tt = dt.*ones(Float64, Nt)
  tt[1] -= dt
  for i = 2:size(t,1)
    tt[t[i,1]] += t[i,2]
  end
  cumsum!(tt, tt)
  return tt
end


text:
# Copyright 2021 Hironori Uchikawa
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Galois Field over 2^8
"""
module GF256

export GF, bvec

const prim_poly = 285 # 1+x^2+x^3+x^4+x^8
const gf_size = 256
const gf_size_m1 = gf_size-1

struct GF
    x::Int64
end

"""
transform a GF symbol to corresponding binary vector
>> bvec(GF(3)) == [1,1,0,0,0,0,0,0]
"""
bvec(a::GF) = digits(Int64, a.x, base=2, pad=8)

Base.:(==)(a::GF, b::GF)::Bool = a.x == b.x

Base.:(+)(a::GF) = a
Base.:(+)(a::GF, b::GF)::GF = GF(a.x ⊻ b.x)

Base.:-(a::GF) = a
Base.:-(a::GF, b::GF)::GF = GF(a.x ⊻ b.x)

Base.:⊻(a::GF) = a
Base.:⊻(a::GF, b::GF)::GF = GF(a.x ⊻ b.x)

"""
makepowertable()

Retrun a list of α^p over GF(256).
"""
function makepowertable()::Array{Int64}
    t = ones(Int64, gf_size)
    v = 1
    for i in 2:gf_size
        v <<= 1
        if v > (gf_size_m1)
            v ⊻= prim_poly
        end
        t[i] = v
    end
    return t
end

"""
powertable[p] is α^p over GF(256)
"""
const powertable = Dict{Int64, Int64}(zip(0:(gf_size_m1), makepowertable()))

"""
logarithm table for GF(256).
logtable[p] means log(p)
"""
const logtable = Dict{Int64, Int64}(zip(makepowertable(), 0:(gf_size-2)))

function Base.log(a::GF)::Int64
    if a.x == 0
        return Inf
    end
    try
        return logtable[a.x]
    catch e
        println("GF(256) symbol must be less than 256!")
        bt = backtrace()
        msg = sprint(showerror, e, bt)
        println(msg)
    end
end

function makemultable()::Array{Int64, 2}
    t = zeros(Int64, gf_size_m1, gf_size_m1)
    for i in 1:gf_size_m1
        for j in 1:gf_size_m1
            t[i,j] = powertable[(logtable[i] + logtable[j]) % (gf_size_m1)]
        end
    end
    return t
end

const multable = makemultable()

function mul(a::GF, b::GF)::GF
    if a.x == 0 | b.x == 0
        return GF(0)
    end
    try
        return GF(multable[a.x,b.x])
    catch e
        println("GF(256) symbol must be less than 256!")
        bt = backtrace()
        msg = sprint(showerror, e, bt)
        println(msg)
    end
end

Base.:*(a::GF, b::GF)::GF = mul(a, b)

function makeinvtable()::Array{Int64}
    t = ones(Int64, gf_size)
    t[1] = 0
    for i in 3:gf_size
        t[i] = powertable[gf_size_m1-logtable[i-1]]
    end
    return t
end

const invtable = Dict{Int64, Int64}(zip(0:(gf_size_m1), makeinvtable()))

function Base.inv(a::GF)::GF
    try
        return GF(invtable[a.x])
    catch e
        println("GF(256) symbol must be less than 256!")
        bt = backtrace()
        msg = sprint(showerror, e, bt)
        println(msg)
    end
end

function makepowtable()::Array{Int64, 2}
    t = ones(Int64, gf_size_m1, gf_size_m1)
    for i in 2:gf_size_m1
        for j in 1:gf_size_m1
            t[i,j] = powertable[(logtable[i] * j) % (gf_size_m1)]
        end
    end
    return t
end

const powtable = makepowtable()

function pow(a::GF, b::Int64)::GF
    if a.x == 0
        return GF(0)
    end
    @assert a.x < gf_size "GF(256) symbol must be less than 256!"
    b %= (gf_size_m1)
    if b == 0
        return GF(1)
    end
    # don't need to think about -b because α^-b is evaluated as inv(α)^b based on literal_pow
    return GF(powtable[a.x,b])
end

Base.:^(a::GF, b::Int64)::GF = pow(a, b)

const α = GF(2)

end # module


text:
using tiddlywinks
using Test

@testset "tiddlywinks.jl" begin
    @test tiddlywinks.tiddly_greet() == "tiddlywinks is working"
    @test tiddlywinks.tiddly_greet() != "hello world"
end


text:
struct FileIntVariate{T<:Signed} <: AbstractIntVariate{T}
    name::String
    length::Int64
    datapath::String
end

Base.length(intvariate::FileIntVariate{T}) where {T<:Signed} = intvariate.length

function slice(intvariate::FileIntVariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Signed}
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    sizeT = sizeof(T)
    iostream = open(intvariate.datapath)
    seek(iostream, sizeT * (fromobs - 1))
    buffer = Vector{T}(undef, slicelength)
    map(Seq(Vector{T}, (buffer, iostream, fromobs, toobs, slicelength), nextdatachunk), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do slice
        view(slice, 1:length(slice))
    end
end

function Base.map(intvariate::FileIntVariate, dataframe::AbstractDataFrame; permute::Bool = false)
    dataframe[getname(intvariate)]
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
export BinaryMetrics

abstract type AbstractBinaryMetric end

struct BinaryMetrics{T <: Real} <: AbstractBinaryMetric
  labels::BitVector
  scores::Vector{T}

  function BinaryMetrics(labels::AbstractVector, scores::AbstractVector)
    L = length(labels)
    scores = float(scores)

    if length(scores) != L
      throw(ErrorException("scores and labels should have the same length."))
    end
    if sum(labels .== true) + sum(labels.==false) != L
      throw(ErrorException("Labels must be either `true` or `false`  (`0` or `1`)"))
    end
    labels = BitArray(labels)
    scores = Array{eltype(scores)}(scores)

    new{eltype(scores)}(labels, scores)
  end
end

"""
  `BinaryMetrics(labels, scores)`

Create an object to compute various metrics of a binary classifier.

For example False Positive Rate at a particular operating point `OP` can be obtained by:
```julia
B = BinaryMetrics(labels, scores)
fpr(B,OP)

```
See documentation for a comprehensive list of available metrics.

"""
BinaryMetrics

Base.length(B::AbstractBinaryMetric) = length(B.labels)
positives(B::AbstractBinaryMetric) = sum(B.labels)
negatives(B::AbstractBinaryMetric) = length(B) - positives(B)

function Base.show(io::IO, B::AbstractBinaryMetric)
  println(io,"Binary Classifier Results ")
  print(io,"Samples:$(length(B)), P:$(positives(B)), N:$(negatives(B))")
end

export positives, negatives
export tp, tn, fp, fn, cm
export tpr, recall, tnr, fpr, fnr
export ppv, fdr
export npv, false_omission_rate, csi
export acc, f1, mcc
export eer, fpr_at_tpr, tpr_at_fpr

is_tp(label, score, OP) =  label && score >= OP
is_tn(label, score, OP) = !label && score  < OP
is_fp(label, score, OP) = !label && score >= OP
is_fn(label, score, OP) =  label && score  < OP

for f in [
          :tp => :is_tp,
          :tn => :is_tn,
          :fp => :is_fp,
          :fn => :is_fn,
         ]
  @eval begin
    function $(f[1])(B::AbstractBinaryMetric, OP::Real)
      cnt = zero(Int)
      for i in eachindex(B.scores)
        if $(f[2])(B.labels[i], B.scores[i], OP)
          cnt += 1
        end
      end
      return cnt
    end

    function $(f[1])(B::AbstractBinaryMetric, OP::AbstractVector)
      cnts = zeros(Int,length(OP))
      for i in eachindex(OP)
        cnts[i] = $(f[1])(B, OP[i])
      end
      return cnts
    end

  end
end

"""
  `tp(B::AbstractBinaryMetric, OP=op(B))`

Returns the true positives of `B` at a particular operating point.
"""
tp

"""
  `tn(B::AbstractBinaryMetric, OP=op(B))`

Returns the true negatives of `B` at a particular operating point.
"""
tn

"""
  `fp(B::AbstractBinaryMetric, OP=op(B))`

Returns the false positives of `B` at a particular operating point.
"""
fp

"""
  `fn(B::AbstractBinaryMetric, OP=op(B))`

Returns the false negatives of `B` at a particular operating point.
"""
fn

"""
  `tpr(B::AbstractBinaryMetric, OP=op(B))`

Returns the true positive rate (aka recall) of `B` at a particular operating point.
"""
tpr(B::AbstractBinaryMetric, OP=op(B)) = tp(B,OP) ./ positives(B)
recall = tpr

"""
  `tnr(B::AbstractBinaryMetric, OP=op(B))`

Returns the true negative rate of `B` at a particular operating point.
"""
tnr(B::AbstractBinaryMetric, OP=op(B)) = tn(B,OP) ./ negatives(B)

"""
  `fnr(B::AbstractBinaryMetric, OP=op(B))`

Returns the false negative rate of `B` at a particular operating point.
"""
fnr(B::AbstractBinaryMetric,OP=op(B)) = fn(B,OP) ./ positives(B)

"""
  `fpr(B::AbstractBinaryMetric, OP=op(B))`

Returns the false positive rate of `B` at a particular operating point.
"""
fpr(B::AbstractBinaryMetric,OP=op(B)) =fp(B,OP) ./ negatives(B)

"""
  `ppv(B::AbstractBinaryMetric, OP=op(B))`

Returns the positive predictive value (aka precision) of `B` at a particular operating point.
"""
function ppv(B::AbstractBinaryMetric,OP=op(B))
  TP = tp(B,OP)
  FP = fp(B,OP)
  return TP ./ ( TP .+ FP )
end

import Base: precision
precision(B::AbstractBinaryMetric,OP=op(B)) = ppv(B,OP)

"""
  `fdr(B::AbstractBinaryMetric, OP=op(B))`

Returns the false discovery rate of `B` at a particular operating point.
"""
function fdr(B::AbstractBinaryMetric,OP=op(B))
  TP = tp(B,OP)
  FP = fp(B,OP)
  return FP ./ ( TP .+ FP )
end

"""
  `npv(B::AbstractBinaryMetric, OP=op(B))`

Returns the negative predictive value of `B` at a particular operating point.
"""
function npv(B::AbstractBinaryMetric, OP=op(B))
  TN = tn(B,OP)
  FN = fn(B,OP)
  return TN ./ ( TN .+ FN )
end

"""
  `false_omission_rate(B::AbstractBinaryMetric, OP=op(B))`

Returns the false omission rate of `B` at a particular operating point.
"""
function false_omission_rate(B::AbstractBinaryMetric, OP=op(B))
  TN = tn(B,OP)
  FN = fn(B,OP)
  return FN ./ ( TN .+ FN )
end

"""
  `csi(B::AbstractBinaryMetric, OP=op(B))`

Returns the critical success index of `B` at a particular operating point.
"""
function csi(B::AbstractBinaryMetric, OP=op(B))
  TP = tp(B,OP)
  FN = fn(B,OP)
  FP = fp(B,OP)
  return TP ./ ( TP .+ FN .+ FP )
end

"""
  `acc(B::AbstractBinaryMetric, OP=op(B))`

Returns the accuracy index of `B` at a particular operating point.
"""
function acc(B::AbstractBinaryMetric, OP=op(B))
  TP = tp(B,OP)
  TN = tn(B,OP)
  return (TP + TN) ./ ( length(B) )
end

"""
  `f1(B::AbstractBinaryMetric, OP=op(B))`

Returns the F1 score of `B` at a particular operating point.
"""
function f1(B::AbstractBinaryMetric, OP=op(B))
  PPV = ppv(B, OP)
  TPR = tpr(B, OP)
  return 2 .*(PPV .* TPR) ./ ( PPV .+ TPR )
end

"""
  `mcc(B::AbstractBinaryMetric, OP=op(B))`

Returns the Matthews correlation coefficient of `B` at a particular operating point.
"""
function mcc(B::AbstractBinaryMetric, OP=op(B))
  TP, FP, FN, TN = tp(B,OP), fp(B,OP), fn(B,OP), tn(B,OP)
  return  @. (TP*TN-FP*FN)/sqrt((TP+FP)*(TP+FN)*(TN+FP)*(TN+FN))
end

"""
  `cm(B::AbstractBinaryMetric, OP::Real)`

Returns the confusion matrix of `B` at a given operating point.
The matrix is defined as follow:

| true positives  | false positives |
| false negatives | true negatives  |
"""
cm(B::AbstractBinaryMetric, OP::Real) =  [ tp(B,OP) fp(B,OP);
                                          fn(B,OP) tn(B,OP) ]

"""
  `cm(B::AbstractBinaryMetric, OP::Array)`

Returns the confusion matrix of `B`. The output is conveniently given as:

`(tps, fps, fns, tns)`

Each element of the Tuple is an array of `length(OP)`.
"""
cm(B::AbstractBinaryMetric, OP::AbstractArray) =  (tp(B,OP), fp(B,OP),
                                                   fn(B,OP), tn(B,OP))
export op

"""
  `op(B::AbstractBinaryMetric)`

Returns significant operating points of `B`.
"""
function op(B::AbstractBinaryMetric)

  OP = sort(unique(B.scores); rev=true) #threshold

  return OP

end

"""
  `eer(FPR::Array,FNR::Array) -> (EER, index)`

Returns the Equal Error Rate (EER) from the intersection of the False Positive rate (must be monotonically increasing) and the False Negative Rate (must be monotonically increasing) and the index of this intersection.

"""
function eer(FPR,FNR)
  if length(FPR) != length(FNR)
    throw(ErrorException("FPR and FNR must have the same length."))
  end
  if !issorted(FPR)
    throw(ErrorException("FPR must be monotonically decreasing."))
  end
  if !issorted(FNR,rev=true)
    throw(ErrorException("FNR must be monotonically increasing."))
  end
  idx = findfirst(FPR .>= FNR)
  EER = (FPR[idx]+FNR[idx])/2
  return EER, idx
end

"""
  `eer(B::AbstractBinaryMetric)`

Returns the Equal Error Rate of `B`.
"""
function eer(B::AbstractBinaryMetric)
  OP = op(B)
  FPR = fpr(B,OP)
  FNR = fnr(B,OP)
  eer(FPR,FNR)
end

"""
  `fpr_at_tpr(B::AbstractBinaryMetric; tpr_at=0.95)`

Returns the false positive rate of `B` when the true positive rate is at `tpr_at` and the corresponding operating point.
"""
function fpr_at_tpr(B::AbstractBinaryMetric, tpr_at=0.95)
  idx = findfirst(tpr(B) .>= tpr_at)
  return fpr(B)[idx], op(B)[idx]
end

"""
  `tpr_at_fpr(B::AbstractBinaryMetric; fpr_at=0.05)`

Returns the true positive rate of `B` when the false positive rate is at `tpr_at` and the corresponding operating point.
"""
function tpr_at_fpr(B::AbstractBinaryMetric, fpr_at=0.05)
  idx = findfirst(fpr(B) .>= fpr_at)
  return tpr(B)[idx], op(B)[idx]
end


text:
"""
    FunctionField{X, Y, Z, C, F, G} <: AbstractField{X, Y, Z, F, G}

An `AbstractField` that returns a function evaluated at location `(X, Y, Z)` (and time, if
`C` is not `Nothing`) when indexed at `i, j, k`.
"""
struct FunctionField{X, Y, Z, C, P, F, G} <: AbstractField{X, Y, Z, F, G}
          func :: F
          grid :: G
         clock :: C
    parameters :: P

    """
        FunctionField{X, Y, Z}(func, grid; clock=nothing, parameters=nothing) where {X, Y, Z}

    Returns a `FunctionField` on `grid` and at location `X, Y, Z`.

    If `clock` is not specified, then `func` must be a function with signature
    `func(x, y, z)`. If clock is specified, `func` must be a function with signature
    `func(x, y, z, t)`, where `t` is internally determined from `clock.time`.

    A FunctionField will return the result of `func(x, y, z [, t])` at `X, Y, Z` on
    `grid` when indexed at `i, j, k`.
    """
    function FunctionField{X, Y, Z}(func, grid; clock=nothing, parameters=nothing) where {X, Y, Z}
        return new{X, Y, Z, typeof(clock),
                   typeof(parameters), typeof(func), typeof(grid)}(func, grid, clock, parameters)
    end

    """
        FunctionField{X, Y, Z}(func::FunctionField, grid; clock) where {X, Y, Z}

    Adds `clock` to an existing `FunctionField` and relocates it to `(X, Y, Z)` on `grid`.
    """
    function FunctionField{X, Y, Z}(f::FunctionField, grid; clock=nothing) where {X, Y, Z}
        return new{X, Y, Z, typeof(clock),
                   typeof(f.parameters), typeof(f.func), typeof(grid)}(f.func, grid, clock, f.parameters)
    end
end

"""Return `a`, or convert `a` to `FunctionField` if `a::Function`"""
fieldify(L, a, grid) = a
fieldify(L, a::Function, grid) = FunctionField(L, a, grid)

"""
    FunctionField(L::Tuple, func, grid)

Returns a stationary `FunctionField` on `grid` and at location `L = (X, Y, Z)`,
where `func` is callable with signature `func(x, y, z)`.
"""
FunctionField(L::Tuple, func, grid) = FunctionField{L[1], L[2], L[3]}(func, grid)

# Ordinary functions needed for fields
architecture(::FunctionField) = nothing
Base.parent(f::FunctionField) = f

# Various possibilities
@inline call_func(clock, parameters, func, x, y, z)     = func(x, y, z, clock.time, parameters)
@inline call_func(::Nothing, parameters, func, x, y, z) = func(x, y, z, parameters)
@inline call_func(clock, ::Nothing, func, x, y, z)      = func(x, y, z, clock.time)
@inline call_func(::Nothing, ::Nothing, func, x, y, z)  = func(x, y, z)

@inline Base.getindex(f::FunctionField{X, Y, Z}, i, j, k) where {X, Y, Z} =
    call_func(f.clock, f.parameters, f.func,
              xnode(X, i, f.grid), ynode(Y, j, f.grid), znode(Z, k, f.grid))

@inline (f::FunctionField)(x, y, z) = call_func(f.clock, f.parameters, f.func, x, y, z)

# set! for function fields
set!(u, f::FunctionField) = set!(u, (x, y, z) -> f.func(x, y, z, f.clock.time, f.parameters))
set!(u, f::FunctionField{X, Y, Z, <:Nothing}) where {X, Y, Z} = set!(u, (x, y, z) -> f.func(x, y, z, f.parameters))
set!(u, f::FunctionField{X, Y, Z, C, <:Nothing}) where {X, Y, Z, C} = set!(u, (x, y, z) -> f.func(x, y, z, f.clock.time))
set!(u, f::FunctionField{X, Y, Z, <:Nothing, <:Nothing}) where {X, Y, Z} = set!(u, (x, y, z) -> f.func(x, y, z))

Adapt.adapt_structure(to, f::FunctionField{X, Y, Z}) where {X, Y, Z} =
    FunctionField{X, Y, Z}(Adapt.adapt(to, f.func),
                           f.grid,
                           clock = Adapt.adapt(to, f.clock),
                           parameters = Adapt.adapt(to, f.parameters))

Base.show(io::IO, field::FunctionField) =
    print(io, "FunctionField located at ", show_location(field), '\n',
          "├── func: $(short_show(field.func))", '\n',
          "├── grid: $(short_show(field.grid))\n",
          "├── clock: $(short_show(field.clock))\n",
          "└── parameters: $(field.parameters)")


text:
function contract!(
  ::Val{:TBLIS},
  R::DenseTensor{ElT},
  labelsR,
  T1::DenseTensor{ElT},
  labelsT1,
  T2::DenseTensor{ElT},
  labelsT2,
  α::ElT,
  β::ElT,
) where {ElT<:LinearAlgebra.BlasReal}
  # TBLIS Tensors
  R_tblis = TBLIS.TTensor{ElT}(array(R), β)
  T1_tblis = TBLIS.TTensor{ElT}(array(T1), α)
  T2_tblis = TBLIS.TTensor{ElT}(array(T2))

  function label_to_char(label)
    # Start at 'a'
    char_start = Char(96)
    if label < 0
      # Start at 'z'
      char_start = Char(123)
    end
    return char_start + label
  end

  function labels_to_tblis(labels)
    if isempty(labels)
      return ""
    end
    str = prod(label_to_char.(labels))
    return str
  end

  labelsT1_tblis = labels_to_tblis(labelsT1)
  labelsT2_tblis = labels_to_tblis(labelsT2)
  labelsR_tblis = labels_to_tblis(labelsR)

  TBLIS.mul!(R_tblis, T1_tblis, T2_tblis, labelsT1_tblis, labelsT2_tblis, labelsR_tblis)

  return R
end


text:
# Sequence rules over the [`SimData`](@ref) object,
# calling [`maprule!`](@ref) for each individual `Rule`.
function sequencerules!(simdata::AbstractSimData)
    newsimdata = sequencerules!(simdata, rules(simdata))
    _maybemask!(grids(newsimdata))
    newsimdata
end
function sequencerules!(simdata::AbstractSimData, rules::Tuple)
    # Mask writes to dest if a mask is provided, except for
    # CellRule which doesn't move values into masked areas
    rule = rules[1]
    rest = tail(rules)
    # Run the first rules
    newsimdata = maprule!(simdata, rule)
    # Run the rest of the rules recursively
    sequencerules!(newsimdata, rest)
end
sequencerules!(simdata::AbstractSimData, rules::Tuple{}) = simdata


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
reject.jl defines methods for the reject(x, y) function
"""

# --- Exports

export reject

# --- Method definitions

using LinearAlgebra: ⋅

"""
    reject(vectors::Matrix, B::AbstractBlade; normalize=false)::Matrix

Compute rejections of `vectors` from `B`. When `normalize` is true, the
rejection vectors are normalized.
"""
function reject end

# ------ Specializations involving a Blade instance

# vectors::Matrix, B::Blade
function reject(vectors::Matrix, B::Blade; normalize::Bool=false)
    # --- Check arguments

    if size(vectors, 1) != dim(B)
        throw(DimensionMismatch("`dim(vectors)` not equal to `dim(B)`"))
    end

    # --- Compute rejections using modified Gram-Schmidt algorithm

    # Initialize rejections
    rejections = Matrix{typeof(volume(B))}(vectors)

    # Remove basis(B) from rejections
    for idx_B in 1:grade(B)
        B_column = basis(B)[:, idx_B]
        for idx in 1:size(rejections, 2)
            rejections[:, idx] -=
                (rejections[:, idx] ⋅ B_column) * B_column
        end
    end

    # Normalize rejection vectors
    if normalize
        for idx in 1:size(rejections, 2)
            norm = LinearAlgebra.norm(rejections[:, idx])
            if norm > 0
                rejections[:, idx] /= norm
            end
        end
    end

    return rejections
end


text:
# Product between PolyVar and Monomial -> Monomial
function Base.:(*)(x::PolyVar{true}, y::PolyVar{true})
    if x === y
        Monomial{true}([x], [2])
    else
        Monomial{true}(x > y ? [x,y] : [y,x], [1,1])
    end
end
function multiplyvar(v::Vector{PolyVar{true}}, x::PolyVar{true}, z)
    i = findfirst(Base.Fix2(<=, x), v)
    if (i != nothing && i > 0) && v[i] == x
        copy(v), multiplyexistingvar(i, z)
    else
        j = (i == nothing || i == 0) ? length(v)+1 : i
        insertvar(v, x, j), insertvar(z, x, j)
    end
end
function Base.:(*)(x::PolyVar{true}, y::Monomial{true})
    w, z = multiplyvar(y.vars, x, y.z)
    Monomial{true}(w, z)
end
Base.:(*)(y::MonomialVector{true}, x::PolyVar{true}) = x * y
function Base.:(*)(x::PolyVar{true}, y::MonomialVector{true})
    w, Z = multiplyvar(y.vars, x, y.Z)
    MonomialVector{true}(w, Z)
end
function _operate_exponents_to!(output::Vector{Int}, op::F, z1::Vector{Int}, z2::Vector{Int}) where {F<:Function}
    @. output = op(z1, z2)
    return
end
function _operate_exponents_to!(output::Vector{Int}, op::F, z1::Vector{Int}, z2::Vector{Int}, maps) where {F<:Function}
    I = maps[1]; i = 1; lI = length(I)
    J = maps[2]; j = 1; lJ = length(J)
    while i <= lI || j <= lJ
        if i > lI || (j <= lJ && J[j] < I[i])
            output[J[j]] = op(0, z2[j])
            j += 1
        elseif j > lJ || (i <= lI && I[i] < J[j])
            output[I[i]] = op(z1[i], 0)
            i += 1
        else
            @assert I[i] == J[j]
            output[I[i]] = op(z1[i], z2[j])
            i += 1
            j += 1
        end
    end
    return
end
# Not used yet
#function _operate_exponents!(op::F, Z::Vector{Vector{Int}}, z2::Vector{Int}, args::Vararg{Any,N}) where {F<:Function,N}
#    return Vector{Int}[_operate_exponents!(op, z, z2, args...) for z in Z]
#end
function _multdivmono!(output, output_variables::Vector{PolyVar{true}},
                      v::Vector{PolyVar{true}}, x::Monomial{true}, op, z)
    if v == x.vars
        if output_variables != v
            resize!(output_variables, length(v))
            copyto!(output_variables, v)
            resize!(output, length(output_variables))
        end
        _operate_exponents_to!(output, op, z, x.z)
    else
        maps = mergevars_to!(output_variables, [v, x.vars])
        resize!(output, length(output_variables))
        _operate_exponents_to!(output, op, z, x.z, maps)
    end
    return
end
function _operate_exponents(op::F, z1::Vector{Int}, z2::Vector{Int}) where {F<:Function}
    return op.(z1, z2)
end
function _operate_exponents(op::F, z1::Vector{Int}, z2::Vector{Int}, n, maps) where {F<:Function}
    newz = zeros(Int, n)
    I = maps[1]; i = 1; lI = length(I)
    J = maps[2]; j = 1; lJ = length(J)
    while i <= lI || j <= lJ
        if i > lI || (j <= lJ && J[j] < I[i])
            newz[J[j]] = op(0, z2[j])
            j += 1
        elseif j > lJ || (i <= lI && I[i] < J[j])
            newz[I[i]] = op(z1[i], 0)
            i += 1
        else
            @assert I[i] == J[j]
            newz[I[i]] = op(z1[i], z2[j])
            i += 1
            j += 1
        end
    end
    return newz
end
function _operate_exponents(op::F, Z::Vector{Vector{Int}}, z2::Vector{Int}, args::Vararg{Any,N}) where {F<:Function,N}
    return Vector{Int}[_operate_exponents(op, z, z2, args...) for z in Z]
end
function multdivmono(v::Vector{PolyVar{true}}, x::Monomial{true}, op, z)
    if v == x.vars
        w = copy(v)
        z_new = _operate_exponents(op, z, x.z)
    else
        w, maps = mergevars([v, x.vars])
        z_new = _operate_exponents(op, z, x.z, length(w), maps)
    end
    return w, z_new
end
function MP.mapexponents_to!(output::Monomial{true}, f::Function, x::Monomial{true}, y::Monomial{true})
    if x.vars == y.vars
        if output.vars != x.vars
            n = length(x.vars)
            resize!(output.vars, n)
            copyto!(output.vars, x.vars)
            resize!(output.z, n)
        end
        _operate_exponents_to!(output.z, f, x.z, y.z)
    else
        _multdivmono!(output.z, output.vars, x.vars, y, f, x.z)
    end
    return output
end
function MP.mapexponents!(f::Function, x::Monomial{true}, y::Monomial{true})
    if x.vars == y.vars
        _operate_exponents_to!(x.z, f, x.z, y.z)
    else
        _multdivmono!(x.z, x.vars, copy(x.vars), y, f, copy(x.z))
    end
    return x
end
function MP.mapexponents(f::Function, x::Monomial{true}, y::Monomial{true})
    w, z = multdivmono(x.vars, y, f, x.z)
    return Monomial{true}(w, z)
end
function Base.:(*)(x::Monomial{true}, y::MonomialVector{true})
    w, Z = multdivmono(y.vars, x, +, y.Z)
    return MonomialVector{true}(w, Z)
end
Base.:(*)(y::MonomialVector{true}, x::Monomial{true}) = x * y
Base.:(*)(x::Monomial{true}, y::PolyVar{true}) = y * x

max_stars_repo_name: blegat/DynamicPolynomials

max_issues_repo_name: blegat/DynamicPolynomials

max_forks_repo_name: blegat/DynamicPolynomials


text:
#=
This file implements functions shared by all discrete spaces.
Discrete spaces are by definition spaces with a finite amount of possible positions.

All these functions are granted "for free" to discrete spaces by simply extending:
- positions(model)
- ids_in_position(position, model)
=#

export positions,
    ids_in_position, agents_in_position, empty_positions, random_empty, has_empty_positions

"""
    positions(model::ABM{<:DiscreteSpace}) → ns
Return an iterator over all positions of a model with a discrete space.

    positions(model::ABM{<:DiscreteSpace}, by::Symbol) → ns
Return all positions of a model with a discrete space, sorting them
using the argument `by` which can be:
* `:random` - randomly sorted
* `:population` - positions are sorted depending on how many agents they accommodate.
  The more populated positions are first.
"""
function positions(model::ABM{<:DiscreteSpace}, by::Symbol)
    n = collect(positions(model))
    itr = reshape(n, length(n))
    if by == :random
        shuffle!(model.rng, itr)
    elseif by == :population
        sort!(itr, by = i -> length(ids_in_position(i, model)), rev = true)
    else
        error("unknown `by`")
    end
    return itr
end

"""
    ids_in_position(position, model::ABM{<:DiscreteSpace})
    ids_in_position(agent, model::ABM{<:DiscreteSpace})

Return the ids of agents in the position corresponding to `position` or position
of `agent`.
"""
ids_in_position(agent::A, model) where {A<:AbstractAgent} =
    ids_in_position(agent.pos, model)

"""
    agents_in_position(position, model::ABM{<:DiscreteSpace})
    agents_in_position(agent, model::ABM{<:DiscreteSpace})

Return the agents in the position corresponding to `position` or position of `agent`.
"""
agents_in_position(agent::A, model) where {A<:AbstractAgent} =
    agents_in_position(agent.pos, model)
agents_in_position(pos, model) = (model[id] for id in ids_in_position(pos, model))

"""
    empty_positions(model)

Return a list of positions that currently have no agents on them.
"""
function empty_positions(model::ABM{<:DiscreteSpace})
    Iterators.filter(i -> length(ids_in_position(i, model)) == 0, positions(model))
end

"""
    isempty(position, model::ABM{<:DiscreteSpace})
Return `true` if there are no agents in `position`.
"""
Base.isempty(pos, model::ABM) = isempty(ids_in_position(pos, model))

"""
    has_empty_positions(model::ABM{<:DiscreteSpace})
Return `true` if there are any positions in the model without agents.
"""
function has_empty_positions(model::ABM{<:DiscreteSpace})
    return any(i -> length(i) == 0, model.space.s)
end

"""
    random_empty(model::ABM{<:DiscreteSpace})
Return a random position without any agents, or `nothing` if no such positions exist.
"""
function random_empty(model::ABM{<:DiscreteSpace})
    empty = empty_positions(model)
    isempty(empty) && return nothing
    rand(model.rng, collect(empty))
end

#######################################################################################
# Discrete space extra agent adding stuff
#######################################################################################
export add_agent_single!, fill_space!, move_agent_single!

"""
    add_agent_single!(agent, model::ABM{<:DiscreteSpace}) → agent

Add the `agent` to a random position in the space while respecting a maximum of one agent
per position, updating the agent's position to the new one.

This function does nothing if there aren't any empty positions.
"""
function add_agent_single!(agent::A, model::ABM{<:DiscreteSpace,A}) where {A<:AbstractAgent}
    position = random_empty(model)
    isnothing(position) && return agent
    agent.pos = position
    add_agent_pos!(agent, model)
    return agent
end

"""
    add_agent_single!(model::ABM{<:DiscreteSpace}, properties...; kwargs...)
Same as `add_agent!(model, properties...; kwargs...)` but ensures that it adds an agent
into a position with no other agents (does nothing if no such position exists).
"""
function add_agent_single!(model::ABM{<:DiscreteSpace}, properties...; kwargs...)
    empty = collect(empty_positions(model))
    if length(empty) > 0
        add_agent!(rand(model.rng, empty), model, properties...; kwargs...)
    end
end

"""
    fill_space!([A ,] model::ABM{<:DiscreteSpace,A}, args...; kwargs...)
    fill_space!([A ,] model::ABM{<:DiscreteSpace,A}, f::Function; kwargs...)
Add one agent to each position in the model's space. Similarly with [`add_agent!`](@ref),
the function creates the necessary agents and
the `args...; kwargs...` are propagated into agent creation.
If instead of `args...` a function `f` is provided, then `args = f(pos)` is the result of
applying `f` where `pos` is each position (tuple for grid, index for graph).

An optional first argument is an agent **type** to be created, and targets mixed agent
models where the agent constructor cannot be deduced (since it is a union).

Example usage in [Daisyworld](@ref).
"""
fill_space!(model::ABM{S,A}, args...; kwargs...) where {S,A<:AbstractAgent} =
    fill_space!(A, model, args...; kwargs...)

function fill_space!(
    ::Type{A},
    model::ABM{<:DiscreteSpace,U},
    args...;
    kwargs...,
) where {A<:AbstractAgent,U<:AbstractAgent}
    for p in positions(model)
        id = nextid(model)
        add_agent_pos!(A(id, p, args...; kwargs...), model)
    end
    return model
end

function fill_space!(
    ::Type{A},
    model::ABM{<:DiscreteSpace,U},
    f::Function;
    kwargs...,
) where {A<:AbstractAgent,U<:AbstractAgent}
    for p in positions(model)
        id = nextid(model)
        args = f(p)
        add_agent_pos!(A(id, p, args...; kwargs...), model)
    end
    return model
end

"""
    move_agent_single!(agent, model::ABM{<:DiscreteSpace}) → agent

Move agent to a random position while respecting a maximum of one agent
per position. If there are no empty positions, the agent won't move.
"""
function move_agent_single!(
    agent::A,
    model::ABM{<:DiscreteSpace,A},
) where {A<:AbstractAgent}
    empty = collect(empty_positions(model))
    if length(empty) > 0
        move_agent!(agent, rand(model.rng, empty), model)
    end
    return agent
end


text:
module CommunalHelperSurfaceSoundPanel

using ..Ahorn, Maple
using Ahorn.CommunalHelper

const directions = String["Up", "Left", "Right"]

@mapdef Entity "CommunalHelper/SurfaceSoundPanel" Panel(
    x::Integer,
    y::Integer,
    soundIndex::Integer=11,
)

const placements = Ahorn.PlacementDict(
    "Sound Surface Panel ($dir) (Communal Helper)" => Ahorn.EntityPlacement(
        Panel,
        "rectangle",
        Dict{String,Any}(
            "orientation" => dir,
        ),
    ) for dir in directions
)

Ahorn.editingOptions(entity::Panel) = Dict{String,Any}(
    "orientation" => Maple.spike_directions,
    "soundIndex" => CommunalHelper.surfaceSounds,
)

Ahorn.minimumSize(entity::Panel) = 8, 8

const resizeDirections = Dict{String,Tuple{Bool,Bool}}(
    "Up" => (true, false),
    "Left" => (false, true),
    "Right" => (false, true),
)

function Ahorn.resizable(entity::Panel)
    orientation = get(entity.data, "orientation", "Up")
    return resizeDirections[orientation]
end

function Ahorn.selection(entity::Panel)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Panel)
    orientation = get(entity.data, "orientation", "Up")

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    x, y = Ahorn.position(entity)
    ox, oy = x + width * (orientation == "Right"), y + height * (orientation == "Down")
    cx, cy = x + width * in(orientation, ("Up", "Down", "Right")),
    y + height * in(orientation, ("Down", "Left", "Right"))

    Ahorn.drawRectangle(ctx, x, y, width, height, (1.0, 0.5, 0.0, 0.4))
    Ahorn.drawLines(ctx, ((ox, oy), (cx, cy)), (1.0, 1.0, 1.0, 1.0), thickness=1)
end

end