text:
struct MyTheme<:InteractBase.WidgetTheme; end
InteractBase.registertheme!(:mytheme, MyTheme())

@testset "theme" begin
    @test gettheme() == NativeHTML()
    settheme!(MyTheme())
    @test gettheme() == MyTheme()
    resettheme!()
    @test gettheme() == NativeHTML()
    settheme!("mytheme")
    @test gettheme() == MyTheme()
    settheme!(:nativehtml)
    @test gettheme() == NativeHTML()
    @test availablethemes() == [:mytheme, :nativehtml]
    @test_throws ErrorException settheme!("not a theme")
end


text:
module DiscFESpacesTests

using Test
using Gridap

D = 2
dom = fill(4,D)
model = CartesianDiscreteModel(partition=tuple(dom...))

order = 1
reffe = PDiscRefFE(Float64,D,order)

fespace = DiscFESpace(reffe,model)

ufun(x) = x[1] + x[2]

trian = Triangulation(model)

nfree = 48
ndiri = 0

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

bh = FEBasis(fespace)
uh = zero(fespace)

cellmat = integrate(inner(bh,bh),trian,quad)
cellvec = integrate(inner(bh,uh),trian,quad)

test_fe_space(fespace, nfree, ndiri, cellmat, cellvec, ufun)

@test celldofids(fespace).data == 1:nfree

uh = interpolate(fespace,ufun)

e = ufun - uh

el2 = sqrt(sum(integrate(inner(e,e),trian,quad)))
@test el2 < 1.0e-10

#ufun(x) = sin(x[1]) * cos(x[2])
#
#uh = interpolate(fespace,ufun)
#
#writevtk(trian,"trian",nref=4,cellfields=["uh"=>uh])

end # module

ext: jl


text:
# test consistency of multithreading option

using Pseudospectra, LinearAlgebra

@testset "Multithreading" begin
@info "Testing multihreading with nt=$(Threads.nthreads())"
for n in [20,60]
    A = randn(n,n) + 1im * randn(n,n)
    Tschur,U,eigA  = schur(A)
    ax = [-3,3,-3,3]
    npts = 100
    opts = Dict{Symbol,Any}()
    Z0,x,y,levels,err,Tproj,eigAproj = psa_compute(Tschur,npts,
                                                   ax,eigA,opts)
    opts = Dict{Symbol,Any}(:threaded => true)
    Z1,x,y,levels,err,Tproj,eigAproj = psa_compute(Tschur,npts,
                                                   ax,eigA,opts)
    @test norm(Z0 - Z1) / norm(Z0) < 1e-4
end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "OpenFOAM"
version = v"8.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/OpenFOAM/OpenFOAM-8/archive/version-8.tar.gz",
                  "94ba11cbaaa12fbb5b356e01758df403ac8832d69da309a5d79f76f42eb008fc"),
    DirectorySource("./bundled"),
]

# In order to set up OpenFOAM, we need to know the version of some of the
# dependencies.
const SCOTCH_VERSION = "6.1.0"

# Bash recipe for building across all platforms
script = "SCOTCH_VERSION=$(SCOTCH_VERSION)\n" * raw"""
cd ${WORKSPACE}/srcdir/OpenFOAM*
atomic_patch -p1 ../patches/etc-bashrc.patch

# Set rpath-link in all C/C++ compilers
LDFLAGS=""
for dir in "" "/dummy" "/mpi-system"; do
    LDFLAGS="${LDFLAGS} -Wl,-rpath-link=${PWD}/platforms/linux64GccDPInt32Opt/lib${dir}"
done
sed -i "s?-m64?-m64 ${LDFLAGS}?g" wmake/rules/*/c*

# Set version of Scotch
echo "export SCOTCH_VERSION=${SCOTCH_VERSION}" > etc/config.sh/scotch
echo "export SCOTCH_ARCH_PATH=${prefix}"      >> etc/config.sh/scotch

# Set up to use our MPI (MPICH)
sed -i 's/WM_MPLIB=SYSTEMOPENMPI/WM_MPLIB=SYSTEMMPI/g' etc/bashrc
export MPI_ROOT="${prefix}"
export MPI_ARCH_FLAGS=""
export MPI_ARCH_INC="-I${includedir}"
export MPI_ARCH_LIBS="-L${libdir} -lmpi"

# Set up the environment.  Note, this script may internally have some failing command, which
# would spring our traps, so we have to allow failures, sigh
source etc/bashrc || true

# Build!
./Allwmake -j${nproc}

# Highly advanced installation process (inspired by Debian:
# https://salsa.debian.org/science-team/openfoam/-/tree/master/debian)
mkdir -p "${libdir}" "${bindir}" "${prefix}/share/openfoam"
cp platforms/linux64GccDPInt32Opt/lib/{,dummy/,mpi-system/}*.${dlext}* "${libdir}/."
cp platforms/linux64GccDPInt32Opt/bin/* "${bindir}/."
cp -r etc/ "${prefix}/share/openfoam/."
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libatmosphericModels", :libatmosphericModels; dont_dlopen=true),
    LibraryProduct("libbarotropicCompressibilityModel", :libbarotropicCompressibilityModel; dont_dlopen=true),
    LibraryProduct("libblockMesh", :libblockMesh; dont_dlopen=true),
    LibraryProduct("libchemistryModel", :libchemistryModel; dont_dlopen=true),
    LibraryProduct("libcoalCombustion", :libcoalCombustion; dont_dlopen=true),
    LibraryProduct("libcombustionModels", :libcombustionModels; dont_dlopen=true),
    LibraryProduct("libconversion", :libconversion; dont_dlopen=true),
    LibraryProduct("libdecompose", :libdecompose; dont_dlopen=true),
    LibraryProduct("libdecompositionMethods", :libdecompositionMethods; dont_dlopen=true),
    LibraryProduct("libdistributed", :libdistributed; dont_dlopen=true),
    LibraryProduct("libdistributionModels", :libdistributionModels; dont_dlopen=true),
    LibraryProduct("libDPMMomentumTransportModels", :libDPMMomentumTransportModels; dont_dlopen=true),
    LibraryProduct("libdriftFluxRelativeVelocityModels", :libdriftFluxRelativeVelocityModels; dont_dlopen=true),
    LibraryProduct("libdriftFluxTransportModels", :libdriftFluxTransportModels; dont_dlopen=true),
    LibraryProduct("libDSMC", :libDSMC; dont_dlopen=true),
    LibraryProduct("libdynamicFvMesh", :libdynamicFvMesh; dont_dlopen=true),
    LibraryProduct("libdynamicMesh", :libdynamicMesh; dont_dlopen=true),
    LibraryProduct("libengine", :libengine; dont_dlopen=true),
    LibraryProduct("libeulerianInterfacialCompositionModels", :libeulerianInterfacialCompositionModels; dont_dlopen=true),
    LibraryProduct("libeulerianInterfacialModels", :libeulerianInterfacialModels; dont_dlopen=true),
    LibraryProduct("libextrude2DMesh", :libextrude2DMesh; dont_dlopen=true),
    LibraryProduct("libextrudeModel", :libextrudeModel; dont_dlopen=true),
    LibraryProduct("libfieldFunctionObjects", :libfieldFunctionObjects; dont_dlopen=true),
    LibraryProduct("libfileFormats", :libfileFormats; dont_dlopen=true),
    LibraryProduct("libfiniteVolume", :libfiniteVolume; dont_dlopen=true),
    LibraryProduct("libfluidThermoMomentumTransportModels", :libfluidThermoMomentumTransportModels; dont_dlopen=true),
    LibraryProduct("libfluidThermophysicalModels", :libfluidThermophysicalModels; dont_dlopen=true),
    LibraryProduct("libfoamToVTK", :libfoamToVTK; dont_dlopen=true),
    LibraryProduct("libforces", :libforces; dont_dlopen=true),
    LibraryProduct("libfvMotionSolvers", :libfvMotionSolvers; dont_dlopen=true),
    LibraryProduct("libfvOptions", :libfvOptions; dont_dlopen=true),
    LibraryProduct("libgenericPatchFields", :libgenericPatchFields; dont_dlopen=true),
    LibraryProduct("libimmiscibleIncompressibleTwoPhaseMixture", :libimmiscibleIncompressibleTwoPhaseMixture; dont_dlopen=true),
    LibraryProduct("libincompressibleMomentumTransportModels", :libincompressibleMomentumTransportModels; dont_dlopen=true),
    LibraryProduct("libincompressibleTransportModels", :libincompressibleTransportModels; dont_dlopen=true),
    LibraryProduct("libincompressibleTwoPhaseMixture", :libincompressibleTwoPhaseMixture; dont_dlopen=true),
    LibraryProduct("libinterfaceProperties", :libinterfaceProperties; dont_dlopen=true),
    LibraryProduct("liblagrangianFunctionObjects", :liblagrangianFunctionObjects; dont_dlopen=true),
    LibraryProduct("liblagrangianIntermediate", :liblagrangianIntermediate; dont_dlopen=true),
    LibraryProduct("liblagrangian", :liblagrangian; dont_dlopen=true),
    LibraryProduct("liblagrangianSpray", :liblagrangianSpray; dont_dlopen=true),
    LibraryProduct("liblagrangianTurbulence", :liblagrangianTurbulence; dont_dlopen=true),
    LibraryProduct("liblaminarFlameSpeedModels", :liblaminarFlameSpeedModels; dont_dlopen=true),
    LibraryProduct("libmeshTools", :libmeshTools; dont_dlopen=true),
    LibraryProduct("libmetisDecomp", :libmetisDecomp; dont_dlopen=true),
    LibraryProduct("libMGridGen", :libMGridGen; dont_dlopen=true),
    LibraryProduct("libmolecularMeasurements", :libmolecularMeasurements; dont_dlopen=true),
    LibraryProduct("libmolecule", :libmolecule; dont_dlopen=true),
    LibraryProduct("libmomentumTransportModels", :libmomentumTransportModels; dont_dlopen=true),
    LibraryProduct("libmultiphaseEulerFoamFunctionObjects", :libmultiphaseEulerFoamFunctionObjects; dont_dlopen=true),
    LibraryProduct("libmultiphaseInterFoam", :libmultiphaseInterFoam; dont_dlopen=true),
    LibraryProduct("libmultiphaseMixtureThermo", :libmultiphaseMixtureThermo; dont_dlopen=true),
    LibraryProduct("libmultiphaseMomentumTransportModels", :libmultiphaseMomentumTransportModels; dont_dlopen=true),
    LibraryProduct("libmultiphaseSystems", :libmultiphaseSystems; dont_dlopen=true),
    LibraryProduct("libmultiphaseThermophysicalTransportModels", :libmultiphaseThermophysicalTransportModels; dont_dlopen=true),
    LibraryProduct("libODE", :libODE; dont_dlopen=true),
    LibraryProduct("libOpenFOAM", :libOpenFOAM; dont_dlopen=true),
    LibraryProduct("libpairPatchAgglomeration", :libpairPatchAgglomeration; dont_dlopen=true),
    LibraryProduct("libphaseChangeTwoPhaseMixtures", :libphaseChangeTwoPhaseMixtures; dont_dlopen=true),
    LibraryProduct("libphaseSystem", :libphaseSystem; dont_dlopen=true),
    LibraryProduct("libpotential", :libpotential; dont_dlopen=true),
    LibraryProduct("libpsiReactionThermophysicalTransportModels", :libpsiReactionThermophysicalTransportModels; dont_dlopen=true),
    LibraryProduct("libPstream", :libPstream; dont_dlopen=true),
    LibraryProduct("libptscotchDecomp", :libptscotchDecomp; dont_dlopen=true),
    LibraryProduct("libradiationModels", :libradiationModels; dont_dlopen=true),
    LibraryProduct("librandomProcesses", :librandomProcesses; dont_dlopen=true),
    LibraryProduct("libreactionThermophysicalModels", :libreactionThermophysicalModels; dont_dlopen=true),
    LibraryProduct("libreconstruct", :libreconstruct; dont_dlopen=true),
    LibraryProduct("libregionModels", :libregionModels; dont_dlopen=true),
    LibraryProduct("librenumberMethods", :librenumberMethods; dont_dlopen=true),
    LibraryProduct("librhoCentralFoam", :librhoCentralFoam; dont_dlopen=true),
    LibraryProduct("librhoReactionThermophysicalTransportModels", :librhoReactionThermophysicalTransportModels; dont_dlopen=true),
    LibraryProduct("librigidBodyDynamics", :librigidBodyDynamics; dont_dlopen=true),
    LibraryProduct("librigidBodyMeshMotion", :librigidBodyMeshMotion; dont_dlopen=true),
    LibraryProduct("librigidBodyState", :librigidBodyState; dont_dlopen=true),
    LibraryProduct("libsampling", :libsampling; dont_dlopen=true),
    LibraryProduct("libscotchDecomp", :libscotchDecomp; dont_dlopen=true),
    LibraryProduct("libsixDoFRigidBodyMotion", :libsixDoFRigidBodyMotion; dont_dlopen=true),
    LibraryProduct("libsixDoFRigidBodyState", :libsixDoFRigidBodyState; dont_dlopen=true),
    LibraryProduct("libSLGThermo", :libSLGThermo; dont_dlopen=true),
    LibraryProduct("libsnappyHexMesh", :libsnappyHexMesh; dont_dlopen=true),
    LibraryProduct("libsolidDisplacementThermo", :libsolidDisplacementThermo; dont_dlopen=true),
    LibraryProduct("libsolidParticle", :libsolidParticle; dont_dlopen=true),
    LibraryProduct("libsolidThermo", :libsolidThermo; dont_dlopen=true),
    LibraryProduct("libsolverFunctionObjects", :libsolverFunctionObjects; dont_dlopen=true),
    LibraryProduct("libspecie", :libspecie; dont_dlopen=true),
    LibraryProduct("libspecieTransfer", :libspecieTransfer; dont_dlopen=true),
    LibraryProduct("libsurfaceFilmDerivedFvPatchFields", :libsurfaceFilmDerivedFvPatchFields; dont_dlopen=true),
    LibraryProduct("libsurfaceFilmModels", :libsurfaceFilmModels; dont_dlopen=true),
    LibraryProduct("libsurfMesh", :libsurfMesh; dont_dlopen=true),
    LibraryProduct("libtabulatedWallFunctions", :libtabulatedWallFunctions; dont_dlopen=true),
    LibraryProduct("libthermalBaffleModels", :libthermalBaffleModels; dont_dlopen=true),
    LibraryProduct("libthermophysicalProperties", :libthermophysicalProperties; dont_dlopen=true),
    LibraryProduct("libthermophysicalTransportModels", :libthermophysicalTransportModels; dont_dlopen=true),
    LibraryProduct("libtopoChangerFvMesh", :libtopoChangerFvMesh; dont_dlopen=true),
    LibraryProduct("libtriSurface", :libtriSurface; dont_dlopen=true),
    LibraryProduct("libtwoPhaseMixture", :libtwoPhaseMixture; dont_dlopen=true),
    LibraryProduct("libtwoPhaseMixtureThermo", :libtwoPhaseMixtureThermo; dont_dlopen=true),
    LibraryProduct("libtwoPhaseProperties", :libtwoPhaseProperties; dont_dlopen=true),
    LibraryProduct("libtwoPhaseSurfaceTension", :libtwoPhaseSurfaceTension; dont_dlopen=true),
    LibraryProduct("libuserd-foam", :libuserd_foam; dont_dlopen=true),
    LibraryProduct("libutilityFunctionObjects", :libutilityFunctionObjects; dont_dlopen=true),
    LibraryProduct("libVoFphaseCompressibleMomentumTransportModels", :libVoFphaseCompressibleMomentumTransportModels; dont_dlopen=true),
    LibraryProduct("libwaves", :libwaves; dont_dlopen=true),
    ExecutableProduct("adiabaticFlameT", :adiabaticFlameT),
    ExecutableProduct("adjointShapeOptimizationFoam", :adjointShapeOptimizationFoam),
    ExecutableProduct("ansysToFoam", :ansysToFoam),
    ExecutableProduct("applyBoundaryLayer", :applyBoundaryLayer),
    ExecutableProduct("attachMesh", :attachMesh),
    ExecutableProduct("autoPatch", :autoPatch),
    ExecutableProduct("autoRefineMesh", :autoRefineMesh),
    ExecutableProduct("blockMesh", :blockMesh),
    ExecutableProduct("boundaryFoam", :boundaryFoam),
    ExecutableProduct("boxTurb", :boxTurb),
    ExecutableProduct("buoyantPimpleFoam", :buoyantPimpleFoam),
    ExecutableProduct("buoyantSimpleFoam", :buoyantSimpleFoam),
    ExecutableProduct("cavitatingFoam", :cavitatingFoam),
    ExecutableProduct("cfx4ToFoam", :cfx4ToFoam),
    ExecutableProduct("changeDictionary", :changeDictionary),
    ExecutableProduct("checkMesh", :checkMesh),
    ExecutableProduct("chemFoam", :chemFoam),
    ExecutableProduct("chemkinToFoam", :chemkinToFoam),
    ExecutableProduct("chtMultiRegionFoam", :chtMultiRegionFoam),
    ExecutableProduct("coalChemistryFoam", :coalChemistryFoam),
    ExecutableProduct("coldEngineFoam", :coldEngineFoam),
    ExecutableProduct("collapseEdges", :collapseEdges),
    ExecutableProduct("combinePatchFaces", :combinePatchFaces),
    ExecutableProduct("compressibleInterFilmFoam", :compressibleInterFilmFoam),
    ExecutableProduct("compressibleInterFoam", :compressibleInterFoam),
    ExecutableProduct("compressibleMultiphaseInterFoam", :compressibleMultiphaseInterFoam),
    ExecutableProduct("createBaffles", :createBaffles),
    ExecutableProduct("createExternalCoupledPatchGeometry", :createExternalCoupledPatchGeometry),
    ExecutableProduct("createPatch", :createPatch),
    ExecutableProduct("datToFoam", :datToFoam),
    ExecutableProduct("decomposePar", :decomposePar),
    ExecutableProduct("deformedGeom", :deformedGeom),
    ExecutableProduct("dnsFoam", :dnsFoam),
    ExecutableProduct("DPMFoam", :DPMFoam),
    ExecutableProduct("driftFluxFoam", :driftFluxFoam),
    ExecutableProduct("dsmcFoam", :dsmcFoam),
    ExecutableProduct("dsmcInitialise", :dsmcInitialise),
    ExecutableProduct("electrostaticFoam", :electrostaticFoam),
    ExecutableProduct("engineCompRatio", :engineCompRatio),
    ExecutableProduct("engineFoam", :engineFoam),
    ExecutableProduct("engineSwirl", :engineSwirl),
    ExecutableProduct("equilibriumCO", :equilibriumCO),
    ExecutableProduct("equilibriumFlameT", :equilibriumFlameT),
    ExecutableProduct("extrude2DMesh", :extrude2DMesh),
    ExecutableProduct("extrudeMesh", :extrudeMesh),
    ExecutableProduct("extrudeToRegionMesh", :extrudeToRegionMesh),
    ExecutableProduct("faceAgglomerate", :faceAgglomerate),
    ExecutableProduct("financialFoam", :financialFoam),
    ExecutableProduct("fireFoam", :fireFoam),
    ExecutableProduct("flattenMesh", :flattenMesh),
    ExecutableProduct("fluent3DMeshToFoam", :fluent3DMeshToFoam),
    ExecutableProduct("fluentMeshToFoam", :fluentMeshToFoam),
    ExecutableProduct("foamDataToFluent", :foamDataToFluent),
    ExecutableProduct("foamDictionary", :foamDictionary),
    ExecutableProduct("foamFormatConvert", :foamFormatConvert),
    ExecutableProduct("foamListTimes", :foamListTimes),
    ExecutableProduct("foamMeshToFluent", :foamMeshToFluent),
    ExecutableProduct("foamSetupCHT", :foamSetupCHT),
    ExecutableProduct("foamToEnsight", :foamToEnsight),
    ExecutableProduct("foamToEnsightParts", :foamToEnsightParts),
    ExecutableProduct("foamToGMV", :foamToGMV),
    ExecutableProduct("foamToStarMesh", :foamToStarMesh),
    ExecutableProduct("foamToSurface", :foamToSurface),
    ExecutableProduct("foamToTetDualMesh", :foamToTetDualMesh),
    ExecutableProduct("foamToVTK", :foamToVTK),
    ExecutableProduct("gambitToFoam", :gambitToFoam),
    ExecutableProduct("gmshToFoam", :gmshToFoam),
    ExecutableProduct("icoFoam", :icoFoam),
    ExecutableProduct("ideasUnvToFoam", :ideasUnvToFoam),
    ExecutableProduct("insideCells", :insideCells),
    ExecutableProduct("interFoam", :interFoam),
    ExecutableProduct("interMixingFoam", :interMixingFoam),
    ExecutableProduct("interPhaseChangeFoam", :interPhaseChangeFoam),
    ExecutableProduct("kivaToFoam", :kivaToFoam),
    ExecutableProduct("laplacianFoam", :laplacianFoam),
    ExecutableProduct("magneticFoam", :magneticFoam),
    ExecutableProduct("mapFields", :mapFields),
    ExecutableProduct("mapFieldsPar", :mapFieldsPar),
    ExecutableProduct("mdEquilibrationFoam", :mdEquilibrationFoam),
    ExecutableProduct("mdFoam", :mdFoam),
    ExecutableProduct("mdInitialise", :mdInitialise),
    ExecutableProduct("mergeMeshes", :mergeMeshes),
    ExecutableProduct("mergeOrSplitBaffles", :mergeOrSplitBaffles),
    ExecutableProduct("mhdFoam", :mhdFoam),
    ExecutableProduct("mirrorMesh", :mirrorMesh),
    ExecutableProduct("mixtureAdiabaticFlameT", :mixtureAdiabaticFlameT),
    ExecutableProduct("modifyMesh", :modifyMesh),
    ExecutableProduct("moveDynamicMesh", :moveDynamicMesh),
    ExecutableProduct("moveEngineMesh", :moveEngineMesh),
    ExecutableProduct("moveMesh", :moveMesh),
    ExecutableProduct("MPPICFoam", :MPPICFoam),
    ExecutableProduct("mshToFoam", :mshToFoam),
    ExecutableProduct("multiphaseEulerFoam", :multiphaseEulerFoam),
    ExecutableProduct("multiphaseInterFoam", :multiphaseInterFoam),
    ExecutableProduct("netgenNeutralToFoam", :netgenNeutralToFoam),
    ExecutableProduct("noise", :noise),
    ExecutableProduct("nonNewtonianIcoFoam", :nonNewtonianIcoFoam),
    ExecutableProduct("objToVTK", :objToVTK),
    ExecutableProduct("orientFaceZone", :orientFaceZone),
    ExecutableProduct("particleFoam", :particleFoam),
    ExecutableProduct("particleTracks", :particleTracks),
    ExecutableProduct("patchSummary", :patchSummary),
    ExecutableProduct("pdfPlot", :pdfPlot),
    ExecutableProduct("PDRFoam", :PDRFoam),
    ExecutableProduct("PDRMesh", :PDRMesh),
    ExecutableProduct("pimpleFoam", :pimpleFoam),
    ExecutableProduct("pisoFoam", :pisoFoam),
    ExecutableProduct("plot3dToFoam", :plot3dToFoam),
    ExecutableProduct("polyDualMesh", :polyDualMesh),
    ExecutableProduct("porousSimpleFoam", :porousSimpleFoam),
    ExecutableProduct("postChannel", :postChannel),
    ExecutableProduct("postProcess", :postProcess),
    ExecutableProduct("potentialFoam", :potentialFoam),
    ExecutableProduct("potentialFreeSurfaceFoam", :potentialFreeSurfaceFoam),
    ExecutableProduct("reactingFoam", :reactingFoam),
    ExecutableProduct("reactingParcelFoam", :reactingParcelFoam),
    ExecutableProduct("reconstructPar", :reconstructPar),
    ExecutableProduct("reconstructParMesh", :reconstructParMesh),
    ExecutableProduct("redistributePar", :redistributePar),
    ExecutableProduct("refineHexMesh", :refineHexMesh),
    ExecutableProduct("refinementLevel", :refinementLevel),
    ExecutableProduct("refineMesh", :refineMesh),
    ExecutableProduct("refineWallLayer", :refineWallLayer),
    ExecutableProduct("removeFaces", :removeFaces),
    ExecutableProduct("renumberMesh", :renumberMesh),
    ExecutableProduct("rhoCentralFoam", :rhoCentralFoam),
    ExecutableProduct("rhoParticleFoam", :rhoParticleFoam),
    ExecutableProduct("rhoPimpleFoam", :rhoPimpleFoam),
    ExecutableProduct("rhoPorousSimpleFoam", :rhoPorousSimpleFoam),
    ExecutableProduct("rhoReactingBuoyantFoam", :rhoReactingBuoyantFoam),
    ExecutableProduct("rhoReactingFoam", :rhoReactingFoam),
    ExecutableProduct("rhoSimpleFoam", :rhoSimpleFoam),
    ExecutableProduct("rotateMesh", :rotateMesh),
    ExecutableProduct("sammToFoam", :sammToFoam),
    ExecutableProduct("scalarTransportFoam", :scalarTransportFoam),
    ExecutableProduct("selectCells", :selectCells),
    ExecutableProduct("setFields", :setFields),
    ExecutableProduct("setSet", :setSet),
    ExecutableProduct("setsToZones", :setsToZones),
    ExecutableProduct("setWaves", :setWaves),
    ExecutableProduct("shallowWaterFoam", :shallowWaterFoam),
    ExecutableProduct("simpleFoam", :simpleFoam),
    ExecutableProduct("simpleReactingParcelFoam", :simpleReactingParcelFoam),
    ExecutableProduct("singleCellMesh", :singleCellMesh),
    ExecutableProduct("smapToFoam", :smapToFoam),
    ExecutableProduct("snappyHexMesh", :snappyHexMesh),
    ExecutableProduct("solidDisplacementFoam", :solidDisplacementFoam),
    ExecutableProduct("solidEquilibriumDisplacementFoam", :solidEquilibriumDisplacementFoam),
    ExecutableProduct("splitCells", :splitCells),
    ExecutableProduct("splitMesh", :splitMesh),
    ExecutableProduct("splitMeshRegions", :splitMeshRegions),
    ExecutableProduct("sprayFoam", :sprayFoam),
    ExecutableProduct("SRFPimpleFoam", :SRFPimpleFoam),
    ExecutableProduct("SRFSimpleFoam", :SRFSimpleFoam),
    ExecutableProduct("star3ToFoam", :star3ToFoam),
    ExecutableProduct("star4ToFoam", :star4ToFoam),
    ExecutableProduct("steadyParticleTracks", :steadyParticleTracks),
    ExecutableProduct("stitchMesh", :stitchMesh),
    ExecutableProduct("subsetMesh", :subsetMesh),
    ExecutableProduct("surfaceAdd", :surfaceAdd),
    ExecutableProduct("surfaceAutoPatch", :surfaceAutoPatch),
    ExecutableProduct("surfaceBooleanFeatures", :surfaceBooleanFeatures),
    ExecutableProduct("surfaceCheck", :surfaceCheck),
    ExecutableProduct("surfaceClean", :surfaceClean),
    ExecutableProduct("surfaceCoarsen", :surfaceCoarsen),
    ExecutableProduct("surfaceConvert", :surfaceConvert),
    ExecutableProduct("surfaceFeatureConvert", :surfaceFeatureConvert),
    ExecutableProduct("surfaceFeatures", :surfaceFeatures),
    ExecutableProduct("surfaceFind", :surfaceFind),
    ExecutableProduct("surfaceHookUp", :surfaceHookUp),
    ExecutableProduct("surfaceInertia", :surfaceInertia),
    ExecutableProduct("surfaceLambdaMuSmooth", :surfaceLambdaMuSmooth),
    ExecutableProduct("surfaceMeshConvert", :surfaceMeshConvert),
    ExecutableProduct("surfaceMeshConvertTesting", :surfaceMeshConvertTesting),
    ExecutableProduct("surfaceMeshExport", :surfaceMeshExport),
    ExecutableProduct("surfaceMeshImport", :surfaceMeshImport),
    ExecutableProduct("surfaceMeshInfo", :surfaceMeshInfo),
    ExecutableProduct("surfaceMeshTriangulate", :surfaceMeshTriangulate),
    ExecutableProduct("surfaceOrient", :surfaceOrient),
    ExecutableProduct("surfacePointMerge", :surfacePointMerge),
    ExecutableProduct("surfaceRedistributePar", :surfaceRedistributePar),
    ExecutableProduct("surfaceRefineRedGreen", :surfaceRefineRedGreen),
    ExecutableProduct("surfaceSplitByPatch", :surfaceSplitByPatch),
    ExecutableProduct("surfaceSplitByTopology", :surfaceSplitByTopology),
    ExecutableProduct("surfaceSplitNonManifolds", :surfaceSplitNonManifolds),
    ExecutableProduct("surfaceSubset", :surfaceSubset),
    ExecutableProduct("surfaceToPatch", :surfaceToPatch),
    ExecutableProduct("surfaceTransformPoints", :surfaceTransformPoints),
    ExecutableProduct("temporalInterpolate", :temporalInterpolate),
    ExecutableProduct("tetgenToFoam", :tetgenToFoam),
    ExecutableProduct("thermoFoam", :thermoFoam),
    ExecutableProduct("topoSet", :topoSet),
    ExecutableProduct("transformPoints", :transformPoints),
    ExecutableProduct("twoLiquidMixingFoam", :twoLiquidMixingFoam),
    ExecutableProduct("viewFactorsGen", :viewFactorsGen),
    ExecutableProduct("vtkUnstructuredToFoam", :vtkUnstructuredToFoam),
    ExecutableProduct("wallFunctionTable", :wallFunctionTable),
    ExecutableProduct("writeMeshObj", :writeMeshObj),
    ExecutableProduct("XiEngineFoam", :XiEngineFoam),
    ExecutableProduct("XiFoam", :XiFoam),
    ExecutableProduct("zipUpMesh", :zipUpMesh),
    FileProduct("share/openfoam/etc", :openfoam_etc),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("MPICH_jll"),
    Dependency("flex_jll"),
    Dependency("SCOTCH_jll"; compat=SCOTCH_VERSION),
    Dependency("PTSCOTCH_jll"),
    Dependency("METIS_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"5", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#-------------------------------------------------
#  generate MMS using the following code
# using SymPy
# x, y = @vars x y
# u = x*y*(1-y)
# v = x*(1-y)

# ux = diff(u, x)
# uy = diff(u, y)
# vx = diff(v, x)
# vy = diff(v, y)

# ε = [
#     ux
#     vy
#     uy + vx
# ]

# H = domain.elements[1].mat[1].H
# σ = H*ε

# f1 = diff(σ[1], x) + diff(σ[3], y)
# f2 = diff(σ[3], x) + diff(σ[2], y)

using NNFEM
using PoreFlow
using PyPlot

m = 40
n = 40
h = 1/m
domain = example_domain(m, n, h)
x, y = domain.nodes[:,1], domain.nodes[:,2]
out = [@. x*y*(1-y)
      @. x*(1-y)]

EBC = zeros(Int64,domain.nnodes, 2)
g = [(@. x*y*(1-y)) (@. x*(1-y))]
FBC = zeros(Int64,domain.nnodes,2)
f = zeros(domain.nnodes,2)

for j = 1:n+1
    EBC[(j-1)*(m+1)+m+1, :] .= -1
    EBC[(j-1)*(m+1)+1, :] .= -1
end

for i = 1:m+1
    EBC[n*(m+1)+i, :] .= -1
    EBC[i, :] .= -1
end

nodes, elements = domain.nodes, domain.elements

domain = Domain(nodes, elements, 2, EBC, g, FBC, f)
H = domain.elements[1].mat[1].H

EBC_func = nothing
FBC_func = nothing

function body_func(x, y, t)
    f1 = @. -1.48148148148148x - 2.46913580246914
    f2 = @. 2.46913580246914 - 4.93827160493827y
    -[f1 f2]
end
globaldata = GlobalData(missing, missing, missing, missing, domain.neqs, EBC_func, FBC_func, body_func)

#-------------------------------------------------
# PoreFlow Solver
@info "Solving using PoreFlow..."
bd = bcnode("all", m, n, h)
K = compute_fem_stiffness_matrix(H, m, n, h)
K, Kbd = fem_impose_Dirichlet_boundary_condition(K, bd, m, n, h)
ff = Kbd * domain.state[[bd; bd.+domain.nnodes]]

F1 = eval_f_on_gauss_pts((x,y)->-1.48148148148148x - 2.46913580246914, m, n, h)
F2 = eval_f_on_gauss_pts((x,y)-> 2.46913580246914 - 4.93827160493827y, m, n, h)
F = compute_fem_source_term(-F1, -F2, m, n, h)

rhs = F - ff
rhs[[bd; bd.+domain.nnodes]] = out[[bd; bd.+domain.nnodes]]
out2 = K\rhs
figure(figsize=(10,10))
title("PoreFlow")
subplot(221)
visualize_scalar_on_scoped_body(out2[1:domain.nnodes], zeros(size(domain.state)...), domain)
subplot(222)
visualize_scalar_on_scoped_body(out[1:domain.nnodes], zeros(size(domain.state)...), domain)

subplot(223)
visualize_scalar_on_scoped_body(out2[domain.nnodes+1:end], zeros(size(domain.state)...), domain)
subplot(224)
visualize_scalar_on_scoped_body(out[domain.nnodes+1:end], zeros(size(domain.state)...), domain)

#-------------------------------------------------
# NNFEM Solver
@info "Solving using NNFEM..."
globaldata, domain = LinearStaticSolver(globaldata, domain)
figure(figsize=(10,10))
title("NNFEM")
subplot(221)
visualize_scalar_on_scoped_body(domain.state[1:domain.nnodes], zeros(size(domain.state)...), domain)
subplot(222)
visualize_scalar_on_scoped_body(out[1:domain.nnodes], zeros(size(domain.state)...), domain)
subplot(223)
visualize_scalar_on_scoped_body(domain.state[domain.nnodes+1:end], zeros(size(domain.state)...), domain)
subplot(224)
visualize_scalar_on_scoped_body(out[domain.nnodes+1:end], zeros(size(domain.state)...), domain)


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d
include("../LF_common/sparse_LF3d.jl")
include("../LF_common/prec_mkl_ilu0.jl")
include("../LF_common/apply_prec_ilu0.jl")
include("build_nabla2_x.jl")
include("build_nabla2_y.jl")
include("build_nabla2_z.jl")

const ⊗ = kron

function test_main( NN::Array{Int64} )

  println("NN = ", NN)

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]
  Npoints = Nx*Ny*Nz

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  nabla2_x = build_nabla2_x(LF)
  nabla2_y = build_nabla2_y(LF)
  nabla2_z = build_nabla2_z(LF)

  # Using Kronecker product
  ∇2 = get_Laplacian3d_kron(LF)

  #
  srand(123)
  v = rand(Npoints)

  # Test multiplication operation
  #@time Lv_1 = nabla2_x*v + nabla2_y*v + nabla2_z*v
  #@time Lv_2 = ∇2*v

  #@printf("Difference: %18.10e\n", sum(abs(Lv_1 - Lv_2))/Npoints)

  # Test the preconditioner

  prec_x = prec_mkl_ilu0( nabla2_x )
  prec_y = prec_mkl_ilu0( nabla2_y )
  prec_z = prec_mkl_ilu0( nabla2_z )

  pvx = apply_prec_ilu0( prec_x, v )
  pvy = apply_prec_ilu0( prec_y, v )
  pvz = apply_prec_ilu0( prec_z, v )
  Pv_1 = pvx + pvy + pvz
  #
  prec_∇2 = prec_mkl_ilu0( ∇2 )
  Pv_2 = apply_prec_ilu0(prec_∇2, v)

  println(sum(v))

  println(sum(pvx))
  println(sum(pvy))
  println(sum(pvz))

  println(sum(Pv_1))
  println(sum(Pv_2))

end

test_main( [3, 3, 3])


text:
include("include.jl")

data = createStructures()

fillInput!(data)
fillOptions!(data)
fillPreprocessor!(data)

if data.options.enumerate
    explicityEnumeration!(data)
    grapfOutputEnum!(data)
else
    singleModel!(data)
    grapfOutput!(data)
    textOutput!(data)
end

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus


text:
mutable struct Queue{T}
    data::Vector{T}
    function (::Type{Queue{T}})() where T
        new{T}(Vector{T}())
    end
end

isempty(q::Queue) = length(q.data) == 0
length(q::Queue) = length(q.data)
peek(q::Queue) = q.data[1]
enqueue!(q::Queue{T}, value::T) where T = push!(q.data, value)
# Enqueue into kth place.
enqueue!(q::Queue{T}, value::T, k::Integer) where T = insert!(q.data, k+1, value)
dequeue!(q::Queue) = popfirst!(q.data)


text:
using Flux: length, reshape, train!, throttle, @epochs
using OperatorLearning, Flux, MAT

device = cpu;

#=
We would like to implement and train a DeepONet that infers the solution
u(x) of the burgers equation on a grid of 1024 points at time one based
on the initial condition a(x) = u(x,0)
=#

# Read the data from MAT file and store it in a dict
# key "a" is the IC
# key "u" is the desired solution at time 1
vars = matread("burgers_data_R10.mat") |> device

# For trial purposes, we might want to train with different resolutions
# So we sample only every n-th element
subsample = 2^3;

# create the x training array, according to our desired grid size
xtrain = vars["a"][1:1000, 1:subsample:end]' |> device;
# create the x test array
xtest = vars["a"][end-99:end, 1:subsample:end]' |> device;

# Create the y training array
ytrain = vars["u"][1:1000, 1:subsample:end] |> device;
# Create the y test array
ytest = vars["u"][end-99:end, 1:subsample:end] |> device;

# The data is missing grid data, so we create it
# `collect` converts data type `range` into an array
grid = collect(range(0, 1, length=1024))' |> device

# Create the DeepONet:
# IC is given on grid of 1024 points, and we solve for a fixed time t in one
# spatial dimension x, making the branch input of size 1024 and trunk size 1
# We choose GeLU activation for both subnets
model = DeepONet((1024,1024,1024),(1,1024,1024),gelu,gelu) |> device

# We use the ADAM optimizer for training
learning_rate = 0.001
opt = ADAM(learning_rate)

# Specify the model parameters
parameters = params(model)

# The loss function
# We can't use the "vanilla" implementation of the mse here since we have
# two distinct inputs to our DeepONet, so we wrap them into a tuple
loss(xtrain,ytrain,sensor) = Flux.Losses.mse(model(xtrain,sensor),ytrain)

# Define a callback function that gives some output during training
evalcb() = @show(loss(xtest,ytest,grid))
# Print the callback only every 5 seconds
throttled_cb = throttle(evalcb, 5)

# Do the training loop
Flux.@epochs 500 train!(loss, parameters, [(xtrain,ytrain,grid)], opt, cb = evalcb)


text:
module SLICOTtools
# Collection of wrappers extracted from SLCOTMath.jl (generated by Ralph Smith)
# based on the SLICOT_jll library (created by the courtesy of Ralph Smith)
using SLICOT_jll
using LinearAlgebra
using LinearAlgebra: BlasInt

export mb03vd!, mb03vy!, mb03bd!, mb03wd!

function chkargsok(ret::BlasInt)
    if ret < 0
        throw(ArgumentError("invalid argument #$(-ret) to SLICOT call"))
    end
end

const BlasBool = BlasInt

"""
    mb03vd!(n::Integer, p::Integer, ilo::Integer, ihi::Integer, A::Array{Float64, 3}, tau::AbstractMatrix{Float64}) -> info::Int64

Reduce a product of `p` real general matrices `A = A_1*A_2*...*A_p`
to upper Hessenberg form, `H = H_1*H_2*...*H_p`, where `H_1` is
upper Hessenberg, and `H_2`, ..., `H_p` are upper triangular, by using
orthogonal similarity transformations on `A`,

        Q_1' * A_1 * Q_2 = H_1,
        Q_2' * A_2 * Q_3 = H_2,
               ...
        Q_p' * A_p * Q_1 = H_p.

The matrices `A_1`, `A_2`, ..., `A_p` are contained in the 3-dimensional array `A`.
The resulting `H_1`, `H_2`, ..., `H_p` and `Q_1`, `Q_2`, ..., `Q_p` overwrite
`A_1`, `A_2`, ..., `A_p` in `A` and the array `tau`.

See the SLICOT documentation of `MB03VD` for details.
"""
function mb03vd!(n::Integer, p::Integer, ilo::Integer, ihi::Integer,
    a::Array{Float64,3}, tau::AbstractMatrix{Float64})

    lda1 = max(1,stride(a,2))
    lda2 = max(1,stride(a,3)÷lda1)
    ldtau = max(1,stride(tau,2))
    info = Ref{BlasInt}()
    dwork = Vector{Float64}(undef, n)

    ccall((:mb03vd_, libslicot), Cvoid, (Ref{BlasInt}, Ref{BlasInt},
            Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
            Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
            Ptr{BlasInt}), n, p, ilo, ihi, a, lda1, lda2, tau,
            ldtau, dwork, info)
    chkargsok(info[])

    return info[]
end

"""
     mb03vy!(n::Integer, p::Integer, ilo::Integer, ihi::Integer, A::Array{Float64, 3}, tau::AbstractMatrix{Float64}) -> info::Int64

Generate the real orthogonal matrices `Q_1`, `Q_2`, ..., `Q_p`,
which are defined as the product of `ihi-ilo` elementary reflectors
of order `n`, as returned in `A_1`, `A_2`, ..., `A_p` by `mb03vd!`:

     Q_j = H_j(ilo) H_j(ilo+1) . . . H_j(ihi-1).

The 3-dimensional arrays `A` and `tau` contains the information on the employed
elementary reflectors. The resulting `Q_1`, `Q_2`, ..., `Q_p` overwrite `A_1`, `A_2`, ..., `A_p`.

See the SLICOT documentation of `MB03VY` for details.

"""
function mb03vy!(n::Integer, p::Integer, ilo::Integer, ihi::Integer,
    a::Array{Float64,3}, tau::AbstractMatrix{Float64})

    lda1 = max(1,stride(a,2))
    lda2 = max(1,stride(a,3)÷lda1)
    ldtau = max(1,stride(tau,2))
    info = Ref{BlasInt}()
    ldwork = BlasInt(-1)
    dwork = Vector{Float64}(undef, 1)

    local jlres
    for iwq in 1:2
        ccall((:mb03vy_, libslicot), Cvoid, (Ref{BlasInt}, Ref{BlasInt},
            Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
            Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
            Ref{BlasInt}, Ptr{BlasInt}), n, p, ilo, ihi, a, lda1,
            lda2, tau, ldtau, dwork, ldwork, info)
        chkargsok(info[])
        if iwq == 1
            ldwork = BlasInt(real(dwork[1]))
            resize!(dwork, ldwork)
        end
    end

    return info[]
end

"""
    mb03vw!(compq::AbstractChar, triu::AbstractChar, qind::AbstractVector{Int64}, k::Integer, n::Integer, h::Integer,
            ilo::Integer, ihi::Integer, s::AbstractVector{Int64}, a::Array{Float64, 3}, q::Array{Float64, 3},
            liwork::Integer, ldwork::Integer) -> info::Int64

Reduce the generalized matrix product

              s[1]           s[2]                 s[k]
      A[:,:,1]     * A[:,:,2]     * ... * A[:,:,k]

to upper Hessenberg-triangular form, where A is N-by-N-by-K and S
is the signature array with values 1 or -1. The H-th matrix of A
is reduced to upper Hessenberg form while the other matrices are
triangularized.

If `compq = 'U'` or `compq = 'I'`, then the orthogonal factors are
computed and stored in the array `Q` so that for `s[i] = 1`,

                    T
        Q[:,:,i](in)   A[:,:,i](in)   Q[:,:,mod(i,k)+1](in)
                                                            T
    =   Q[:,:,i](out)  A[:,:,i](out)  Q[:,:,mod(i,k)+1](out),

and for `s[i] = -1`,

                             T
        Q[:,:,mod(i,k)+1](in)   A[:,:,i](in)   Q[:,:,i](in)
                                                            T
    =   Q[:,:,mod(i,k)+1](out)  A[:,:,i](out)  Q[:,:,i](out).

A partial generation of the orthogonal factors can be realized
via the array `qind`.

If `triu = 'N'` only matrices with negative signature are reduced to upper
triangular form in the first stage of the algorithm.
If `triu = 'A'` all possible `n-1` matrices with negative signature are reduced.

See the SLICOT documentation of `MB03VW` for details.
"""
function mb03vw!(compq::AbstractChar, qind::AbstractVector{BlasInt}, atriu::AbstractChar,
    n::Integer, k::Integer, h::Integer, ilo::Integer, ihi::Integer,
    s::AbstractVector{BlasInt}, a::Array{Float64,3},
    q::Array{Float64,3}, liwork::Integer, ldwork::Integer)

    lda1 = max(1,stride(a,2))
    lda2 = max(1,stride(a,3)÷lda1)
    ldq1 = max(1,stride(q,2))
    ldq2 = max(1,stride(q,3)÷ldq1)
    info = Ref{BlasInt}()
    iwork = Vector{BlasInt}(undef, liwork)
    dwork = Vector{Float64}(undef, ldwork)
    hi = (h < 0 || h > k) ? 1 : h

    ccall((:mb03vw_, libslicot), Cvoid, (Ref{UInt8}, Ptr{BlasInt},
            Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
            Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{Float64},
            Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
            Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
            Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}, Clong, Clong),
            compq, qind, atriu, n, k, hi,
            ilo, ihi, s, a, lda1, lda2, q, ldq1, ldq2,
            iwork, liwork, dwork, ldwork, info, 1, 1)
    chkargsok(info[])

    return info[]
end

"""
    mb03bd!(job::AbstractChar, defl::AbstractChar, compq::AbstractChar, qind::AbstractVector{Int64}, k::Integer, n::Integer, h::Integer,
            ilo::Integer, ihi::Integer, s::AbstractVector{Int64}, a::Array{Float64, 3}, q::Array{Float64, 3}, alphar::AbstractVector{Float64},
            alphai::AbstractVector{Float64}, beta::AbstractVector{Float64}, scal::AbstractVector{Int64}, liwork::Integer, ldwork::Integer) -> (info::Int64, iwarn::Int64)

Find the eigenvalues of the generalized matrix product

              s[1]           s[2]                 s[k]
      A[:,:,1]     * A[:,:,2]     * ... * A[:,:,k]

where `A[:,:,h]` is upper Hessenberg and `A[:,:,i]`, `i <> h`, is upper
triangular, using a double-shift version of the periodic
QZ method. In addition, `A` may be reduced to periodic Schur form:
`A[:,:,h]` is upper quasi-triangular and all the other factors
`A[:,:,i]` are upper triangular. Optionally, the 2-by-2 triangular
matrices corresponding to 2-by-2 diagonal blocks in `A[:,:,h]`
are so reduced that their product is a 2-by-2 diagonal matrix.

If `compq = 'U'` or `compq = 'I'`, then the orthogonal factors are
computed and stored in the array `Q` so that for `s[i] = 1`,

                    T
        Q[:,:,i](in)   A[:,:,i](in)   Q[:,:,mod(i,k)+1](in)
                                                            T
    =   Q[:,:,i](out)  A[:,:,i](out)  Q[:,:,mod(i,k)+1](out),

and for `s[i] = -1`,

                             T
        Q[:,:,mod(i,k)+1](in)   A[:,:,i](in)   Q[:,:,i](in)
                                                            T
    =   Q[:,:,mod(i,k)+1](out)  A[:,:,i](out)  Q[:,:,i](out).

A partial generation of the orthogonal factors can be realized
via the array `qind`.

See the SLICOT documentation of `MB03BD` for details.
"""
function mb03bd!(job::AbstractChar, defl::AbstractChar,
    compq::AbstractChar, qind::AbstractVector{BlasInt}, k::Integer,
    n::Integer, h::Integer, ilo::Integer, ihi::Integer,
    s::AbstractVector{BlasInt}, a::Array{Float64,3},
    q::Array{Float64,3},
    alphar::AbstractVector{Float64}, alphai::AbstractVector{Float64},
    beta::AbstractVector{Float64}, scal::AbstractVector{BlasInt},
    liwork::Integer, ldwork::Integer)

    lda1 = max(1,stride(a,2))
    lda2 = max(1,stride(a,3)÷lda1)
    ldq1 = max(1,stride(q,2))
    ldq2 = max(1,stride(q,3)÷ldq1)
    info = Ref{BlasInt}()
    iwarn = Ref{BlasInt}()
    iwork = Vector{BlasInt}(undef, liwork)
    dwork = Vector{Float64}(undef, ldwork)

    ccall((:mb03bd_, libslicot), Cvoid, (Ref{UInt8}, Ref{UInt8},
            Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
            Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt},
            Ptr{Float64}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64},
            Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64},
            Ptr{Float64}, Ptr{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
            Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
            Clong, Clong, Clong), job, defl, compq, qind, k, n, h,
            ilo, ihi, s, a, lda1, lda2, q, ldq1, ldq2, alphar,
            alphai, beta, scal, iwork, liwork, dwork, ldwork, iwarn,
            info, 1, 1, 1)
    chkargsok(info[])

    return info[], iwarn[]
end

function mb03bz!(job::AbstractChar, compq::AbstractChar, k::Integer,
    n::Integer, ilo::Integer, ihi::Integer,
    s::AbstractVector{BlasInt}, a::Array{ComplexF64,3}, q::Array{ComplexF64,3},
    alpha::AbstractVector{ComplexF64},
    beta::AbstractVector{ComplexF64}, scal::AbstractVector{BlasInt},
    ldwork::Integer, lzwork::Integer)
"""
    mb03bz!(job::AbstractChar, compq::AbstractChar, k::Integer, n::Integer,
            ilo::Integer, ihi::Integer, s::AbstractVector{Int64}, a::Array{ComplexF64, 3},
            q::Array{ComplexF64, 3}, alpha::AbstractVector{ComplexF64}, beta::AbstractVector{ComplexF64},
            scal::AbstractVector{Int64}, ldwork::Integer, lzwork::Integer) -> info::Int64

"""

    lda1 = max(1,stride(a,2))
    lda2 = max(1,stride(a,3)÷lda1)
    ldq1 = max(1,stride(q,2))
    ldq2 = max(1,stride(q,3)÷ldq1)
    info = Ref{BlasInt}()
    dwork = Vector{Float64}(undef, ldwork)
    zwork = Vector{ComplexF64}(undef, lzwork)

    ccall((:mb03bz_, libslicot), Cvoid, (Ref{UInt8}, Ref{UInt8},
            Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
            Ptr{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt},
            Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt},
            Ref{BlasInt}, Ptr{ComplexF64}, Ptr{ComplexF64},
            Ptr{BlasInt}, Ptr{Float64}, Ref{BlasInt},
            Ptr{ComplexF64}, Ref{BlasInt}, Ptr{BlasInt}, Clong,
            Clong), job, compq, k, n, ilo, ihi, s, a, lda1, lda2, q,
            ldq1, ldq2, alpha, beta, scal, dwork, ldwork, zwork,
            lzwork, info, 1, 1)
    chkargsok(info[])

    return info[]
end
"""
    mb03wd!(job::AbstractChar, compz::AbstractChar, n::Integer, p::Integer,
            ilo::Integer, ihi::Integer, iloz::Integer, ihiz::Integer,  h::Array{Float64, 3}, z::Array{Float64, 3},
            wr::AbstractVector{Float64}, wi::AbstractVector{Float64}, ldwork::Integer) -> info::Int64

Compute the Schur decomposition and the eigenvalues of a
product of matrices, H = `H_1`*`H_2`*...*`H_p`, with `H_1` an upper
Hessenberg matrix and `H_2`, ..., `H_p` upper triangular matrices,
without evaluating the product. Specifically, the matrices Z_i
are computed, such that

        `Z_1' * H_1 * Z_2 = T_1,`
        `Z_2' * H_2 * Z_3 = T_2,`
               `...`
        `Z_p' * H_p * Z_1 = T_p,`

where `T_1` is in real Schur form, and `T_2`, ..., `T_p` are upper
triangular.

The routine works primarily with the Hessenberg and triangular
submatrices in rows and columns ILO to IHI, but optionally applies
the transformations to all the rows and columns of the matrices
H_i, i = 1,...,p. The transformations can be optionally
accumulated.

See the SLICOT documentation of `MB03WD` for details.
"""
function mb03wd!(job::AbstractChar, compz::AbstractChar, n::Integer,
    p::Integer, ilo::Integer, ihi::Integer, iloz::Integer,
    ihiz::Integer, h::Array{Float64,3}, z::Array{Float64,3},
    wr::AbstractVector{Float64}, wi::AbstractVector{Float64},
    ldwork::Integer)

    ldh1 = max(1,stride(h,2))
    ldh2 = max(1,stride(h,3)÷ldh1)
    ldz1 = max(1,stride(z,2))
    ldz2 = max(1,stride(z,3)÷ldz1)
    info = Ref{BlasInt}()
    dwork = Vector{Float64}(undef, ldwork)

    ccall((:mb03wd_, libslicot), Cvoid, (Ref{UInt8}, Ref{UInt8},
            Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
            Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
            Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ref{BlasInt},
            Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{BlasInt},
            Ptr{BlasInt}, Clong, Clong), job, compz, n, p, ilo, ihi,
            iloz, ihiz, h, ldh1, ldh2, z, ldz1, ldz2, wr, wi, dwork,
            ldwork, info, 1, 1)
    chkargsok(info[])

    return info[]
end

end # module


text:
using Nemo, GibbsTypePriors, JLD

grid_k = [collect(1:25:10000); 10000]

Vnk_numerical_accuracy = GibbsTypePriors.Vnk_NGG.(10000, grid_k, 1.2, 0.6) |> x -> accuracy_bits.(x)
save("test/graphical_tests/saves_for_graphical_tests/accuracy_Vnk_10000.jld", "Vnk_numerical_accuracy", Vnk_numerical_accuracy)

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


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

function SCIPhasCurrentNodeLP(scip)
    ccall((:SCIPhasCurrentNodeLP, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPisLPConstructed(scip)
    ccall((:SCIPisLPConstructed, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPconstructLP(scip, cutoff)
    ccall((:SCIPconstructLP, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{UInt32}), scip, cutoff)
end

function SCIPflushLP(scip)
    ccall((:SCIPflushLP, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPgetLPSolstat(scip)
    ccall((:SCIPgetLPSolstat, libscip), SCIP_LPSOLSTAT, (Ptr{SCIP_},), scip)
end

function SCIPisLPPrimalReliable(scip)
    ccall((:SCIPisLPPrimalReliable, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPisLPDualReliable(scip)
    ccall((:SCIPisLPDualReliable, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPisLPRelax(scip)
    ccall((:SCIPisLPRelax, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPgetLPObjval(scip)
    ccall((:SCIPgetLPObjval, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLPColumnObjval(scip)
    ccall((:SCIPgetLPColumnObjval, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLPLooseObjval(scip)
    ccall((:SCIPgetLPLooseObjval, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetGlobalPseudoObjval(scip)
    ccall((:SCIPgetGlobalPseudoObjval, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetPseudoObjval(scip)
    ccall((:SCIPgetPseudoObjval, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPisRootLPRelax(scip)
    ccall((:SCIPisRootLPRelax, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPgetLPRootObjval(scip)
    ccall((:SCIPgetLPRootObjval, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLPRootColumnObjval(scip)
    ccall((:SCIPgetLPRootColumnObjval, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLPRootLooseObjval(scip)
    ccall((:SCIPgetLPRootLooseObjval, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLPColsData(scip, cols, ncols)
    ccall((:SCIPgetLPColsData, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{Ptr{SCIP_COL}}}, Ptr{Cint}), scip, cols, ncols)
end

function SCIPgetLPCols(scip)
    ccall((:SCIPgetLPCols, libscip), Ptr{Ptr{SCIP_COL}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNLPCols(scip)
    ccall((:SCIPgetNLPCols, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetLPRowsData(scip, rows, nrows)
    ccall((:SCIPgetLPRowsData, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{Ptr{SCIP_ROW}}}, Ptr{Cint}), scip, rows, nrows)
end

function SCIPgetLPRows(scip)
    ccall((:SCIPgetLPRows, libscip), Ptr{Ptr{SCIP_ROW}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNLPRows(scip)
    ccall((:SCIPgetNLPRows, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPallColsInLP(scip)
    ccall((:SCIPallColsInLP, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPisLPSolBasic(scip)
    ccall((:SCIPisLPSolBasic, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPgetLPBasisInd(scip, basisind)
    ccall((:SCIPgetLPBasisInd, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cint}), scip, basisind)
end

function SCIPgetLPBInvRow(scip, r, coefs, inds, ninds)
    ccall((:SCIPgetLPBInvRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), scip, r, coefs, inds, ninds)
end

function SCIPgetLPBInvCol(scip, c, coefs, inds, ninds)
    ccall((:SCIPgetLPBInvCol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), scip, c, coefs, inds, ninds)
end

function SCIPgetLPBInvARow(scip, r, binvrow, coefs, inds, ninds)
    ccall((:SCIPgetLPBInvARow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), scip, r, binvrow, coefs, inds, ninds)
end

function SCIPgetLPBInvACol(scip, c, coefs, inds, ninds)
    ccall((:SCIPgetLPBInvACol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), scip, c, coefs, inds, ninds)
end

function SCIPsumLPRows(scip, weights, sumcoef, sumlhs, sumrhs)
    ccall((:SCIPsumLPRows, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cdouble}, Ptr{SCIP_REALARRAY}, Ptr{Cdouble}, Ptr{Cdouble}), scip, weights, sumcoef, sumlhs, sumrhs)
end

function SCIPwriteLP(scip, filename)
    ccall((:SCIPwriteLP, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, filename)
end

function SCIPwriteMIP(scip, filename, genericnames, origobj, lazyconss)
    ccall((:SCIPwriteMIP, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, UInt32, UInt32, UInt32), scip, filename, genericnames, origobj, lazyconss)
end

function SCIPgetLPI(scip, lpi)
    ccall((:SCIPgetLPI, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_LPI}}), scip, lpi)
end

function SCIPprintLPSolutionQuality(scip, file)
    ccall((:SCIPprintLPSolutionQuality, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPcomputeLPRelIntPoint(scip, relaxrows, inclobjcutoff, timelimit, iterlimit, point)
    ccall((:SCIPcomputeLPRelIntPoint, libscip), SCIP_RETCODE, (Ptr{SCIP_}, UInt32, UInt32, Cdouble, Cint, Ptr{Ptr{SCIP_SOL}}), scip, relaxrows, inclobjcutoff, timelimit, iterlimit, point)
end

function SCIPgetColRedcost(scip, col)
    ccall((:SCIPgetColRedcost, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_COL}), scip, col)
end

function SCIPgetColFarkasCoef(scip, col)
    ccall((:SCIPgetColFarkasCoef, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_COL}), scip, col)
end

function SCIPmarkColNotRemovableLocal(scip, col)
    ccall((:SCIPmarkColNotRemovableLocal, libscip), Cvoid, (Ptr{SCIP_}, Ptr{SCIP_COL}), scip, col)
end

function SCIPcreateRowCons(scip, row, conshdlr, name, len, cols, vals, lhs, rhs, _local, modifiable, removable)
    ccall((:SCIPcreateRowCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}, Ptr{SCIP_CONSHDLR}, Cstring, Cint, Ptr{Ptr{SCIP_COL}}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, UInt32, UInt32), scip, row, conshdlr, name, len, cols, vals, lhs, rhs, _local, modifiable, removable)
end

function SCIPcreateRowSepa(scip, row, sepa, name, len, cols, vals, lhs, rhs, _local, modifiable, removable)
    ccall((:SCIPcreateRowSepa, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}, Ptr{SCIP_SEPA}, Cstring, Cint, Ptr{Ptr{SCIP_COL}}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, UInt32, UInt32), scip, row, sepa, name, len, cols, vals, lhs, rhs, _local, modifiable, removable)
end

function SCIPcreateRowUnspec(scip, row, name, len, cols, vals, lhs, rhs, _local, modifiable, removable)
    ccall((:SCIPcreateRowUnspec, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}, Cstring, Cint, Ptr{Ptr{SCIP_COL}}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, UInt32, UInt32), scip, row, name, len, cols, vals, lhs, rhs, _local, modifiable, removable)
end

function SCIPcreateRow(scip, row, name, len, cols, vals, lhs, rhs, _local, modifiable, removable)
    ccall((:SCIPcreateRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}, Cstring, Cint, Ptr{Ptr{SCIP_COL}}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, UInt32, UInt32), scip, row, name, len, cols, vals, lhs, rhs, _local, modifiable, removable)
end

function SCIPcreateEmptyRowCons(scip, row, conshdlr, name, lhs, rhs, _local, modifiable, removable)
    ccall((:SCIPcreateEmptyRowCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}, Ptr{SCIP_CONSHDLR}, Cstring, Cdouble, Cdouble, UInt32, UInt32, UInt32), scip, row, conshdlr, name, lhs, rhs, _local, modifiable, removable)
end

function SCIPcreateEmptyRowSepa(scip, row, sepa, name, lhs, rhs, _local, modifiable, removable)
    ccall((:SCIPcreateEmptyRowSepa, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}, Ptr{SCIP_SEPA}, Cstring, Cdouble, Cdouble, UInt32, UInt32, UInt32), scip, row, sepa, name, lhs, rhs, _local, modifiable, removable)
end

function SCIPcreateEmptyRowUnspec(scip, row, name, lhs, rhs, _local, modifiable, removable)
    ccall((:SCIPcreateEmptyRowUnspec, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}, Cstring, Cdouble, Cdouble, UInt32, UInt32, UInt32), scip, row, name, lhs, rhs, _local, modifiable, removable)
end

function SCIPcreateEmptyRow(scip, row, name, lhs, rhs, _local, modifiable, removable)
    ccall((:SCIPcreateEmptyRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}, Cstring, Cdouble, Cdouble, UInt32, UInt32, UInt32), scip, row, name, lhs, rhs, _local, modifiable, removable)
end

function SCIPcaptureRow(scip, row)
    ccall((:SCIPcaptureRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPreleaseRow(scip, row)
    ccall((:SCIPreleaseRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_ROW}}), scip, row)
end

function SCIPchgRowLhs(scip, row, lhs)
    ccall((:SCIPchgRowLhs, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Cdouble), scip, row, lhs)
end

function SCIPchgRowRhs(scip, row, rhs)
    ccall((:SCIPchgRowRhs, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Cdouble), scip, row, rhs)
end

function SCIPcacheRowExtensions(scip, row)
    ccall((:SCIPcacheRowExtensions, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPflushRowExtensions(scip, row)
    ccall((:SCIPflushRowExtensions, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPaddVarToRow(scip, row, var, val)
    ccall((:SCIPaddVarToRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Ptr{SCIP_VAR}, Cdouble), scip, row, var, val)
end

function SCIPaddVarsToRow(scip, row, nvars, vars, vals)
    ccall((:SCIPaddVarsToRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}), scip, row, nvars, vars, vals)
end

function SCIPaddVarsToRowSameCoef(scip, row, nvars, vars, val)
    ccall((:SCIPaddVarsToRowSameCoef, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Cint, Ptr{Ptr{SCIP_VAR}}, Cdouble), scip, row, nvars, vars, val)
end

function SCIPcalcRowIntegralScalar(scip, row, mindelta, maxdelta, maxdnom, maxscale, usecontvars, intscalar, success)
    ccall((:SCIPcalcRowIntegralScalar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Cdouble, Cdouble, Clonglong, Cdouble, UInt32, Ptr{Cdouble}, Ptr{UInt32}), scip, row, mindelta, maxdelta, maxdnom, maxscale, usecontvars, intscalar, success)
end

function SCIPmakeRowIntegral(scip, row, mindelta, maxdelta, maxdnom, maxscale, usecontvars, success)
    ccall((:SCIPmakeRowIntegral, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Cdouble, Cdouble, Clonglong, Cdouble, UInt32, Ptr{UInt32}), scip, row, mindelta, maxdelta, maxdnom, maxscale, usecontvars, success)
end

function SCIPmarkRowNotRemovableLocal(scip, row)
    ccall((:SCIPmarkRowNotRemovableLocal, libscip), Cvoid, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowNumIntCols(scip, row)
    ccall((:SCIPgetRowNumIntCols, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowMinCoef(scip, row)
    ccall((:SCIPgetRowMinCoef, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowMaxCoef(scip, row)
    ccall((:SCIPgetRowMaxCoef, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowMinActivity(scip, row)
    ccall((:SCIPgetRowMinActivity, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowMaxActivity(scip, row)
    ccall((:SCIPgetRowMaxActivity, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPrecalcRowLPActivity(scip, row)
    ccall((:SCIPrecalcRowLPActivity, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowLPActivity(scip, row)
    ccall((:SCIPgetRowLPActivity, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowLPFeasibility(scip, row)
    ccall((:SCIPgetRowLPFeasibility, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPrecalcRowPseudoActivity(scip, row)
    ccall((:SCIPrecalcRowPseudoActivity, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowPseudoActivity(scip, row)
    ccall((:SCIPgetRowPseudoActivity, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowPseudoFeasibility(scip, row)
    ccall((:SCIPgetRowPseudoFeasibility, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPrecalcRowActivity(scip, row)
    ccall((:SCIPrecalcRowActivity, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowActivity(scip, row)
    ccall((:SCIPgetRowActivity, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowFeasibility(scip, row)
    ccall((:SCIPgetRowFeasibility, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetRowSolActivity(scip, row, sol)
    ccall((:SCIPgetRowSolActivity, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Ptr{SCIP_SOL}), scip, row, sol)
end

function SCIPgetRowSolFeasibility(scip, row, sol)
    ccall((:SCIPgetRowSolFeasibility, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Ptr{SCIP_SOL}), scip, row, sol)
end

function SCIPprintRow(scip, row, file)
    ccall((:SCIPprintRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Ptr{FILE}), scip, row, file)
end

function SCIPstartDive(scip)
    ccall((:SCIPstartDive, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPendDive(scip)
    ccall((:SCIPendDive, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPchgCutoffboundDive(scip, newcutoffbound)
    ccall((:SCIPchgCutoffboundDive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble), scip, newcutoffbound)
end

function SCIPchgVarObjDive(scip, var, newobj)
    ccall((:SCIPchgVarObjDive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, newobj)
end

function SCIPchgVarLbDive(scip, var, newbound)
    ccall((:SCIPchgVarLbDive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, newbound)
end

function SCIPchgVarUbDive(scip, var, newbound)
    ccall((:SCIPchgVarUbDive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, newbound)
end

function SCIPaddRowDive(scip, row)
    ccall((:SCIPaddRowDive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPchgRowLhsDive(scip, row, newlhs)
    ccall((:SCIPchgRowLhsDive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Cdouble), scip, row, newlhs)
end

function SCIPchgRowRhsDive(scip, row, newrhs)
    ccall((:SCIPchgRowRhsDive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, Cdouble), scip, row, newrhs)
end

function SCIPgetVarObjDive(scip, var)
    ccall((:SCIPgetVarObjDive, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarLbDive(scip, var)
    ccall((:SCIPgetVarLbDive, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarUbDive(scip, var)
    ccall((:SCIPgetVarUbDive, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPsolveDiveLP(scip, itlim, lperror, cutoff)
    ccall((:SCIPsolveDiveLP, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{UInt32}, Ptr{UInt32}), scip, itlim, lperror, cutoff)
end

function SCIPgetLastDivenode(scip)
    ccall((:SCIPgetLastDivenode, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPinDive(scip)
    ccall((:SCIPinDive, libscip), UInt32, (Ptr{SCIP_},), scip)
end


text:
using DataStructures
using LightGraphs, SimpleWeightedGraphs

function prim_heuristic(graph::SimpleWeightedGraph,
    k::Int;
    startnode::Union{Nothing,Int}=nothing,
    upperbound=typemax(Int))::Union{Tuple{Vector{Edge{Int}}, Int}, Nothing}
    n = nv(graph)
    es = edges(graph)
    es = collect(SimpleWeightedEdge, es)
    sort!(es, by=weight)
    k -= 1
    if isnothing(startnode)
        candidateedges = []
        candidatenodes = Dict{Int,Int}()
        for e in es
            i = src(e)
            j = dst(e)
            if i > 1 && j > 1
                if isempty(candidateedges)
                    push!(candidateedges, e)
                    push!(candidatenodes, i => 0)
                    push!(candidatenodes, j => 0)
                else
                    w = weight(e)
                    wh = weight(candidateedges[1])
                    if w < wh
                        candidateedges = [e]
                        candidatenodes = Dict{Int,Int}(i => 0, j => 0)
                    elseif w == wh
                        push!(candidateedges, e)
                        push!(candidatenodes, i => 0)
                        push!(candidatenodes, j => 0)
                    end
                end
            end
        end
        @assert !isempty(candidatenodes) "No candidate nodes."
        for e in es
            i = src(e)
            j = dst(e)
            ks = keys(candidatenodes)
            if i in ks || j in ks
                if i in ks
                    candidatenodes[i] = get(candidatenodes, i, 0) + Int(round(weight(e)))
                end
                if j in ks
                    candidatenodes[j] = get(candidatenodes, j, 0) + Int(round(weight(e)))
                end
            end
        end
        (_, startnode) = findmin(candidatenodes)
    end
    @assert !isnothing(startnode) "No startnode"
    @assert startnode != 1 "Startnode is 1"
    @assert 1 <= startnode && startnode <= n "Startnode $startnode is out of bounds (1, $n)"
    time = @elapsed begin
        finished::Vector{Bool} = zeros(Bool, n)
        finished[startnode] = true
        parents::Vector{Int} = zeros(Int, n)

        treesize::Int = 0
        treeweight::Int = 0
        lastnode::Int = -1

        while treesize < k && treeweight < upperbound && !isempty(es)
        del = 0
        for e in es
            del += 1
            u = src(e)
            v = dst(e)
            ((u == 1 || v == 1) || finished[u] == finished[v]) && continue
            w = Int(round(weight(e)))
            if finished[u]
                finished[v] = true
                parents[v] = u
                lastnode = u
            else
                finished[u] = true
                parents[u] = v
                lastnode = v
            end
            treeweight += w
            treesize += 1
            break
        end
        deleteat!(es, del)
    end
    end

    if treeweight > upperbound
        @debug "Prim $startnode: >$upperbound after $treesize nodes, in $(format_seconds_readable(time))."
        return nothing
    end

    kmst::Vector{Edge{Int}} = [Edge{Int}(parents[v], v) for v in vertices(graph) if parents[v] != 0]

    @debug "Prim $startnode: $treeweight, in $(format_seconds_readable(time))."

    return (kmst, treeweight)

end

function kruskal_heuristic(graph::SimpleWeightedGraph, k::Int; unskippable::Union{Nothing,Edge{Int}} = nothing, bound::Int = typemax(Int))::Tuple{Union{Vector{Edge{Int}}, Nothing}, Int, Set{Edge{Int}}}
    time = @elapsed begin
    skips::Int = 0
    skipped :: Set{Edge{Int}} = Set()
    n::Int= nv(graph)
    es::Vector{SimpleWeightedEdge} = collect(SimpleWeightedEdge, edges(graph))
    sort!(es, by=weight)

    intree::Vector{Bool} = zeros(Bool, n)
    parents::Vector{Int} = zeros(Int, n)

    s::Int = 0
    d::Int = 0
    w::Int = 0
    treeweight::Int = 0
    treesize::Int = 0

    if !isnothing(unskippable)
        s = src(unskippable)
        d = dst(unskippable)
        @assert 1 <= s <= n "s: $s is not in [1,n]."
        @assert 1 <= d <= n "d: $d is not in [1,n]."
        intree[s] = true
        intree[d] = true
        parents[d] = s
        treesize += 1
        for e in es
            (src(e) == s && dst(e) == d) || (src(e) == d && dst(e) == s) || continue
            treeweight += Int(round(Int, weight(e)))
        end
    end

    for e::SimpleWeightedEdge in Iterators.drop(es, n-1)
        treesize < (k-1) || break

        s = src(e)
        d = dst(e)

        isnothing(unskippable) || !((src(e) == src(unskippable) && dst(e) == dst(unskippable)) || ((src(e) == dst(unskippable) && dst(e) == src(unskippable)))) || continue

        if intree[s] && intree[d]
            skips += 1
            push!(skipped, Edge{Int}(min(s,d), max(s,d)))
            continue
        end

        if intree[s] || !intree[d]
            parents[d] = s
        else
            parents[s] = d
        end

        treesize += 1

        intree[s] = true
        intree[d] = true
        treeweight += Int(round(Int,weight(e)))
        if treeweight > bound
            break
        end
    end
    end # time-measuring end

    if treeweight > bound
        @debug "Kruskal $(src(unskippable))-$(dst(unskippable)): >$bound after $treesize nodes, in $(format_seconds_readable(time)) with $(skips == 0 ? "no skips" : (skips == 1 ? "1 skip" : "$skips skips"))."
        return (nothing, treeweight, skipped)
    end

    @assert treesize == k-1 "Tree should have $(k-1) edges but has $treesize."

    @debug "Kruskal$(isnothing(unskippable) ? "" : " $(src(unskippable))-$(dst(unskippable))"): $treeweight, in $(format_seconds_readable(time)) with $(skips == 0 ? "no skips" : (skips == 1 ? "1 skip" : "$skips skips"))."

    return ([Edge{Int}(v, parents[v]) for v in vertices(graph) if parents[v] != 0], treeweight, skipped)

end

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
"""
Main code for running the branching process. This invokes other Julia scripts in order
to run:
    - read_in_TP.jl:
        Code for reading the TP csv files saved from the fitting and
        forecasting parts of the model.
    - helper_functions.jl:
        Code which is commonly reused in the model. This involves non-specialised functions
        for sampling onset and infection times, and for sampling from a NegativeBinomial and
        Binomial distribution using asymptotic limiting properties.
    - forecast_types.jl
        Code for some user defined types. This lets us group variables of importance. There
        is a struct for the dynamical constants, properties of the data and a
        simple mapping for moving between types and indices.
    - simulation_consistency.jl
        Code for checking consistency of simulations against the data. The code in this file
        is also responsible for injecting cases into the simulation if we are below some
        threshold.
    - assumptions.jl
        One of the most important files. Keeps track of all assumptions in a function that is
        used to initialise the system constants.

The key logic of this simulation is featured in the function simulate_branching_process(). This
is written in such a way that it is relatively straightforward to follow the logic of the processes
and functions are called in a systematic way. simulate_branching_process() will return reshaped
arrays that are the "good" simulated realisations based off the consistency checks.
"""

using Distributions
using Random
using Statistics
using LinearAlgebra
using ProgressBars
using DataStructures

include("read_in_TP.jl")
include("helper_functions.jl")
include("forecast_types.jl")
include("simulation_consistency.jl")
include("assumptions.jl")

##### Methods

function initialise_state_arrays(
    sim_duration,
    observation_period,
    nsims
)
    """
    Initialising a state array object with arrays of zeros.
    We pad the infection array Z with 35 days to account for
    infections occuring prior to the simulation period. We do this
    separately (and not in a struct) as the arrays are large.
    """

    Z = zeros(Int, sim_duration+35, 3, nsims)
    D = zeros(Int, sim_duration, 3, nsims)
    U = zeros(Int, sim_duration, 3, nsims)

    sim_realisations = Realisations(Z, D, U)

    return sim_realisations
end

function initialise_population!(
    forecast::Forecast,
    start_day,
    D0,
)
    """
    Fill the infection array for each simulation, inititialising the states of
    all the simulations.
    """
    Z = forecast.sim_realisations.Z
    D = forecast.sim_realisations.D
    U = forecast.sim_realisations.U

    # extract prob of being symptomatic and prob of symptoms given detection
    p_symp = forecast.sim_constants.p_symp
    p_detect_given_symp = forecast.sim_constants.p_detect_given_symp
    T_end = forecast.sim_features.T_end
    individual_type_map = forecast.individual_type_map

    # get total initial cases
    total_initial_cases = sum(D0)

    for sim in 1:size(Z, 3)

        # infer some initial undetected symptomatic
        num_symptomatic_undetected = 0
        if D0.S == 0
            num_symptomatic_undetected = sample_negative_binomial_limit(1, p_detect_given_symp)
        else
            num_symptomatic_undetected = sample_negative_binomial_limit(D0.S, p_detect_given_symp)
        end

        # infer some initial undetected asymptomatic
        total_symptomatic = D0.S + num_symptomatic_undetected
        num_asymptomatic_undetected = 0
        if total_symptomatic == 0
            num_asymptomatic_undetected = sample_negative_binomial_limit(1, p_symp)
        else
            num_asymptomatic_undetected = sample_negative_binomial_limit(total_symptomatic, p_symp)
        end

        # total initial cases include the undetected
        total_initial_cases_sim = total_initial_cases + num_symptomatic_undetected + num_asymptomatic_undetected

        if total_initial_cases_sim > 0

            infection_times = zeros(Int, total_initial_cases_sim)

            # sample negative infection times as these individuals are assumed to occur
            # before the simulation
            for i in 1:total_initial_cases_sim
                infection_times[i] = -sample_onset_time()
            end

            # save the initial day upon which to start simulating from
            start_day[sim] = minimum(infection_times)

            # add the individuals to the state arrays
            infection_time_counter = 1

            for _ in 1:D0.S
                Z[infection_times[infection_time_counter]+36,individual_type_map.S,sim] += 1
                infection_time_counter += 1
                D[1,individual_type_map.S,sim] += 1
            end

            for _ in 1:D0.A
                Z[infection_times[infection_time_counter]+36,individual_type_map.A,sim] += 1
                infection_time_counter += 1
                D[1,individual_type_map.A,sim] += 1
            end

            for _ in 1:num_symptomatic_undetected
                Z[infection_times[infection_time_counter]+36,individual_type_map.S,sim] += 1
                infection_time_counter += 1
                U[1,individual_type_map.S,sim] += 1
            end

            for _ in 1:num_asymptomatic_undetected
                Z[infection_times[infection_time_counter]+36,individual_type_map.A,sim] += 1
                infection_time_counter += 1
                U[1,individual_type_map.A,sim] += 1
            end
        end
    end

    return nothing
end

function import_cases_model!(
	forecast::Forecast,
    import_cases,
	forecast_start_date,
)
	"""
	A model to handle international imports in the forecasting of cases. This can be run prior
	to the standard forecasting model and will pre-fill D, U and Z for each simulation.

	We assume that cases
	arise from
	D[t] ∼ NegBin(a[t], 1/(b+1))
	where
	a[t] = α + f(I[t])
	b = β + 1
	are the parameters. The function f() is an exponential weighted moving average such that
	f(I[t]) = ϕ I[t] + (1-ϕ) f(I[t-1])
	and α, β are hyperparameters fixed at 0.5 and 0.2, respectively. We take b = β + 1 as we
	assume a fixed period of 1 day for estimating the posterior.
	"""

    Z = forecast.sim_realisations.Z
    D = forecast.sim_realisations.D
    U = forecast.sim_realisations.U

    # get the number of observation days
    T_observed = forecast.sim_features.T_observed
    # duration of forecast
    T_end = forecast.sim_features.T_end

	# number of days after t = 0 (2020/03/01)
	days_into_forecast = Day(Date(forecast_start_date) - Date("2020-03-01")).value

	# the number of days the simulation is run for
	forecast_days_plus_padding = -10:T_end-5
	# posterior a, b for sampling number of imports on day t
	a_post = zeros(Float64, length(forecast_days_plus_padding))

	# the current ema is taken 11 days before the forecast start date
	current_ema = zero(Float64)
    current_ema += import_cases[days_into_forecast - 11]

	for (i, t) in enumerate(forecast_days_plus_padding)
		# if within observation period, update the posterior. Otherwise use the last
		# available measurement
		if t <= T_observed
			count_on_day = import_cases[days_into_forecast + t]
			current_ema = forecast.sim_constants.ϕ*count_on_day + (1-forecast.sim_constants.ϕ)*current_ema
			a_post[i] = forecast.sim_constants.prior_alpha + current_ema
		else
			a_post[i] = forecast.sim_constants.prior_alpha + current_ema
		end
	end

    # assuming a fixed period of 1 day (so this is the same for the whole period)
    b_post = forecast.sim_constants.prior_beta + 1

    # preallocate daily import vectors
    D_I = zeros(Int, length(a_post))
    U_I = zeros(Int, length(a_post))
    unobserved_I = zeros(Int, length(a_post))

	for sim in 1:size(D, 3)

		# sample the number of imported cases detected and undetected (broadcasting to allocate in place)
		D_I .= rand.(NegativeBinomial.(a_post, b_post / (1 + b_post)))
        for (i, d) in enumerate(D_I)
    		unobserved_I[i] = d == 0 ? 1 : d
        end
		U_I .= rand.(NegativeBinomial.(unobserved_I, forecast.sim_constants.qi))

		for (i, t) in enumerate(forecast_days_plus_padding)
            # day of infection is time t which is why the infection time in
			if D_I[i] + U_I[i] > 0
				# add sampled imports to observation arrays
                assign_to_arrays_and_times!(
                    forecast,
                    t,
                    sim,
                    num_imports_detected=D_I[i],
					num_imports_undetected=U_I[i],
                )
			end
		end
	end

	return nothing

end

function get_proportion_infected(
    day,
    cases_pre_forecast,
    D_sim,
    U_sim,
    N
)
    """
    Calculate the proportion of cases that have been infected up to current day.
    """
    total_cumulative_cases = cases_pre_forecast

    for j in 1:3
        for i in 1:day
            total_cumulative_cases += D_sim[i,j] + U_sim[i,j]
        end
    end

    proportion_infected = total_cumulative_cases / N

    return proportion_infected
end

function sample_offspring!(
    forecast::Forecast,
    t,
    TP_local_sim,
    TP_import_sim,
    TP_indices,
    susceptible_depletion_sim,
    sim,
)
    """
    Sample offspring for all adults at time t. This uses the fact that the sum of N NB(s, p) is
    NB(N*s, p) and hence we can sample ALL offspring for a particular type of adults.
    """
    Z = forecast.sim_realisations.Z
    D = forecast.sim_realisations.D
    U = forecast.sim_realisations.U

    # get the α's
    α_s = forecast.sim_constants.α_s
    α_a = forecast.sim_constants.α_a
    # assumptions surrouding the probability of symptomatic, relative infectiousness γ and the ratio of Reff (α's) (CHECK THIS)
    p_symp = forecast.sim_constants.p_symp
    p_detect_given_symp = forecast.sim_constants.p_detect_given_symp
    p_detect_given_asymp = forecast.sim_constants.p_detect_given_asymp
    k = forecast.sim_constants.k

    # extract fixed simulation values
    cases_pre_forecast = forecast.sim_features.cases_pre_forecast
    N = forecast.sim_features.N

    # initialise number of offspring outside of the loop
    num_offspring = zero(Int)

    # pointer to the current day's infections
    Z_tmp = @view Z[t+36,:,sim]

    # if the total number of infected individuals at time t is 0, return immediately
    if sum(Z_tmp) > 0
        # extract the number of each parent who are infected on day t
        (S_parents, A_parents, I_parents) = Z_tmp

        # reset infections on day t (so we can inject new cases)
        # find day based on the indices representing time from forecast origin
        TP_ind = findfirst(t == ind for ind in TP_indices)

        # create views and use a function barrier to decrease allocations
        D_sim = @views D[:,:,sim]
        U_sim = @views U[:,:,sim]

        # take max of 0 and TP to deal with cases of depletion factor > TP accouting for susceptible_depletion
        proportion_infected = get_proportion_infected(t, cases_pre_forecast, D_sim, U_sim, N)
        TP_local_parent = max(0, TP_local_sim[TP_ind] * (1 - susceptible_depletion_sim*proportion_infected))
        TP_import_parent = max(0, TP_import_sim[TP_ind] * (1 - susceptible_depletion_sim*proportion_infected))

        # total number of infections arising from all parents infected at time t
        # (sum of M NB(r,p) is NB(r*M,p))
        if S_parents > 0 && TP_local_parent > 0
            s = k*S_parents
            p = 1 - α_s*TP_local_parent/(k+α_s*TP_local_parent)
            num_offspring += sample_negative_binomial_limit(s, p)
        end

        if A_parents > 0 && TP_local_parent > 0
            s = k*A_parents
            p = 1 - α_a*TP_local_parent/(k+α_a*TP_local_parent)
            num_offspring += sample_negative_binomial_limit(s, p)
        end

        if I_parents > 0 && TP_import_parent > 0
            s = k*I_parents
            # p_{v,h} is the proportion of hotel quarantine workers vaccinated
            p_vh = 0.9+rand(Beta(2, 4))*9/100
            # v_{e,h} is the overall vaccine effectiveness
            v_eh = 0.83+rand(Beta(2, 2))*14/100
            TP_import_parent *= (1-p_vh*v_eh)*1.39*1.3
            p = 1 - TP_import_parent/(k+TP_import_parent)
            num_offspring += sample_negative_binomial_limit(s, p)
        end

        # we dont ever get import offspring
        if num_offspring > 0
            # sample those symptomatic
            num_symptomatic_offspring = sample_binomial_limit(num_offspring, p_symp)
            num_asymptomatic_offspring = num_offspring-num_symptomatic_offspring
            # sample number detected of each type
            num_symptomatic_offspring_detected = sample_binomial_limit(num_symptomatic_offspring, p_detect_given_symp)
            num_asymptomatic_offspring_detected = sample_binomial_limit(num_asymptomatic_offspring, p_detect_given_asymp)
            # calc number undetected
            num_symptomatic_offspring_undetected = num_symptomatic_offspring-num_symptomatic_offspring_detected
            num_asymptomatic_offspring_undetected = num_asymptomatic_offspring-num_asymptomatic_offspring_detected

            # add the onsets to observation arrays
            assign_to_arrays_and_times!(
                forecast,
                t,
                sim,
                num_symptomatic_detected=num_symptomatic_offspring_detected,
                num_symptomatic_undetected=num_symptomatic_offspring_undetected,
                num_asymptomatic_detected=num_asymptomatic_offspring_detected,
                num_asymptomatic_undetected=num_asymptomatic_offspring_undetected,
            )

            # zero out infections on day t so we can restart sim
            # and keep the progress
            Z_tmp .= 0
        end
    end

    return nothing
end

function assign_to_arrays_and_times!(
    forecast::Forecast,
    day,
    sim;
    num_symptomatic_detected=0,
    num_symptomatic_undetected=0,
    num_asymptomatic_detected=0,
    num_asymptomatic_undetected=0,
    num_imports_detected=0,
    num_imports_undetected=0,
)
    """
    Increment the arrays D and U based on sampled offsrping counts. The number of each offspring
    are keyword parameters which are default 0. This allows us to call this function from multiple
    spots.
    """

    # TODO: refactor this using multiple dispatch.

    Z = forecast.sim_realisations.Z
    D = forecast.sim_realisations.D
    U = forecast.sim_realisations.U
    individual_type_map = forecast.individual_type_map
    omicron_dominant_day = forecast.sim_features.omicron_dominant_day
    T_end = forecast.sim_features.T_end

    # sampling detected cases
    for _ in 1:num_symptomatic_detected
        infection_time = day + sample_infection_time(omicron = day >= omicron_dominant_day)
        if infection_time <= T_end
            Z[infection_time+36,individual_type_map.S,sim] += 1
        end

        onset_time = infection_time + sample_onset_time(omicron = infection_time >= omicron_dominant_day)
        if onset_time < 1
            D[1,individual_type_map.S,sim] += 1
        elseif onset_time <= T_end
            # time+1 as first index is day 0
            D[onset_time+1,individual_type_map.S,sim] += 1
        end
    end

    for _ in 1:num_asymptomatic_detected
        infection_time = day + sample_infection_time(omicron = day >= omicron_dominant_day)
        if infection_time <= T_end
            Z[infection_time+36,individual_type_map.A,sim] += 1
        end

        onset_time = infection_time + sample_onset_time(omicron = infection_time >= omicron_dominant_day)
        if onset_time < 1
            D[1,individual_type_map.A,sim] += 1
        elseif onset_time <= T_end
            # time+1 as first index is day 0
            D[onset_time+1,individual_type_map.A,sim] += 1
        end
    end

    # sampling undetected cases
    for _ in 1:num_symptomatic_undetected
        infection_time = day + sample_infection_time(omicron = day >= omicron_dominant_day)
        if infection_time <= T_end
            Z[infection_time+36,individual_type_map.S,sim] += 1
        end

        onset_time = infection_time + sample_onset_time(omicron = infection_time >= omicron_dominant_day)
        if onset_time < 1
            U[1,individual_type_map.S,sim] += 1
        elseif onset_time <= T_end
            # time+1 as first index is day 0
            U[onset_time+1,individual_type_map.S,sim] += 1
        end
    end

    for _ in 1:num_asymptomatic_undetected
        infection_time = day + sample_infection_time(omicron = day >= omicron_dominant_day)
        if infection_time <= T_end
            Z[infection_time+36,individual_type_map.A,sim] += 1
        end

        onset_time = infection_time + sample_onset_time(omicron = infection_time >= omicron_dominant_day)
        if onset_time < 1
            U[1,individual_type_map.A,sim] += 1
        elseif onset_time <= T_end
            # time+1 as first index is day 0
            U[onset_time+1,individual_type_map.A,sim] += 1
        end
    end

    # imports are treated slightly differently as the infection time is inferred
    # through the simpler model
    for _ in 1:num_imports_detected
        infection_time = day
        if infection_time <= T_end
            Z[infection_time+36,individual_type_map.I,sim] += 1
        end

        onset_time = infection_time + sample_onset_time(omicron = infection_time >= omicron_dominant_day)
        if onset_time < 1
            D[1,individual_type_map.I,sim] += 1
        elseif onset_time <= T_end
            # time+1 as first index is day 0
            D[onset_time+1,individual_type_map.I,sim] += 1
        end
    end

    for _ in 1:num_imports_undetected
        infection_time = day
        if infection_time <= T_end
            Z[infection_time+35,individual_type_map.I,sim] += 1
        end

        onset_time = infection_time + sample_onset_time(omicron = infection_time >= omicron_dominant_day)
        if onset_time < 1
            U[1,individual_type_map.I,sim] += 1
        elseif onset_time <= T_end
            # time+1 as first index is day 0
            U[onset_time+1,individual_type_map.I,sim] += 1
        end
    end

    return nothing
end

function get_simulation_limits(
    local_cases,
    forecast_start_date,
    omicron_dominant_date,
    cases_pre_forecast,
    TP_indices,
    N,
    T_observed,
    T_end
)
    """
    Using the observed cases, determine the limits of cases over the backcast and
    nowcast. This assumes consistency over windows of fixed length and a nowcast period
    of 14 days.
    """

    days_delta = (Dates.Date(omicron_dominant_date) - Dates.Date(forecast_start_date)).value

    cases_pre_backcast = sum(@view local_cases[1:days_delta])
    cases_backcast = sum(@view local_cases[days_delta+1:T_observed])
    cases_pre_nowcast = sum(@view local_cases[21:T_observed])
    cases_nowcast = sum(@view local_cases[10:T_observed])

    # # take the cumulative sum of the local cases for easiness of
    # # calculating the incidence over periods of time
    # cumulative_local_cases = cumsum(local_cases)

    # # extract the backcast and nowcast cases
    # backcast_cases = cumulative_local_cases[1:end-14]
    # backcast_cases = backcast_cases[end:-20:1]
    # nowcast_cases = cumulative_local_cases[end-13:end]
    # nowcast_cases = nowcast_cases[end]
    # # vector of observed cases in each window
    # cumulative_cases_for_windows = [reverse!(backcast_cases); nowcast_cases]
    # cases_in_each_window = Int.(diff([0; cumulative_cases_for_windows]))
    # # fill vector of days each window holds
    # window_lengths = ones(Int,length(cumulative_cases_for_windows))     # NOTE: not used yet

    # for (i, val) in enumerate(cumulative_cases_for_windows)
    #     window_lengths[i] = findlast(cumulative_local_cases .== val)
    # end

    cases_in_each_window = [cases_pre_backcast, cases_backcast, cases_pre_nowcast, cases_nowcast]

    # calculate minimum and maximum observed cases in each period
    min_cases = floor.(Int, [0.3*cases_pre_backcast, 0.5*cases_backcast, 0.5*cases_pre_nowcast, 0.6*cases_nowcast])
    # min_cases = 0*cases_in_each_window
    max_cases = ceil.(Int, [2.5*cases_pre_backcast, 2.5*cases_backcast, 3.0*cases_pre_nowcast, 4.0*cases_nowcast])

    # assume maximum of 250 cases if the observed is less than that
    for (i, val) in enumerate(max_cases)
        if val < 100
            max_cases[i] = 100
        end
    end

    for (i, val) in enumerate(min_cases)
        if val < 50
            min_cases[i] = 0
        end
    end

    # the maximum allowable cases over the forecast period is the population size
    max_forecast_cases = N

    # get the day we want to start using omicron GI and incubation period (+1) as 0 corresponds
    # to the first element of the arrays
    omicron_dominant_day = (omicron_dominant_date - forecast_start_date).value + 1

    sim_features = Features(
        max_forecast_cases,
        cases_pre_forecast,
        N,
        T_observed,
        T_end,
        omicron_dominant_day,
    )

    window_lengths = 0

    return (
        sim_features,
        window_lengths,
        min_cases,
        max_cases,
    )
end

function simulate_branching_process(
    D0,
    N,
    nsims,
    local_cases,
    import_cases,
    cases_pre_forecast,
    forecast_start_date,
    date,
    omicron_dominant_date,
    state,
)
    """
    Simulate branching process nsims times conditional on the cumulative observed cases, with
    initial state X0 = (S, A, I) (a named tuple). T is the duration of the simulation and N is the
    population size which is currently unused.
    """
    # read in the TP for a given state and date
    (TP_indices, TP_local, TP_import) = create_state_TP_matrices(forecast_start_date, date, state)
    # read in susceptible_depletion parameters
    susceptible_depletion = read_in_susceptible_depletion(date)

    # length of observed time series
    T_observed = length(local_cases)
    # duration of forecast simulation
    T_end = sum(ind > 0 for ind in TP_indices)

    cumulative_local_cases = cumsum(local_cases)
    max_restarts = 10
    good_sims = ones(Bool, nsims)
    good_TPs = zeros(Bool, size(TP_local, 2))

    # get simulation constants
    (sim_constants, individual_type_map) = set_simulation_constants(state)
    # calculate the upper and lower limits for checking sims against data
    (sim_features, window_lengths, min_cases, max_cases) = get_simulation_limits(
        local_cases,
        forecast_start_date,
        omicron_dominant_date,
        cases_pre_forecast,
        TP_indices,
        N,
        T_observed,
        T_end,
    )

    # the range of days for the simulation
    day_range = -35:T_end

    # initialise state arrays (T_end + 1 as we include 0 day)
    sim_realisations = initialise_state_arrays(T_end+1, T_observed+1, nsims)

    # put everything into a forecast object
    forecast = Forecast(
        sim_features,
        sim_realisations,
        sim_constants,
        individual_type_map,
    )

    # vector to hold the first day of infections for each simulation.
    start_day = zeros(Int, nsims)

    # generate initial infected population
    initialise_population!(
        forecast,
        start_day,
        D0,
    )

    import_cases_model!(
        forecast,
        import_cases,
        forecast_start_date,
    )

    # for sim in 1:nsims
    # for sim in ProgressBar(1:nsims)
    Threads.@threads for sim in ProgressBar(1:nsims)
        # counts in each window
        case_counts = zeros(Int, length(max_cases))
        # sample the TP/susceptible_depletion for this sim
        TP_ind = sim % 2000 == 0 ? 2000 : sim % 2000
        TP_local_sim = @view TP_local[:,TP_ind]
        TP_import_sim = @view TP_import[:,TP_ind]
        susceptible_depletion_sim = susceptible_depletion[TP_ind]

        # indices for looping through the sorted indices
        bad_sim = false
        injected_cases = false

        # initialise counter to first day
        day = start_day[sim]
        last_injection = 0
        # counter for the number of times we've had to restart the simulation
        n_restarts = 0

        while day <= day_range[end]

            # simulate offspring of parents
            sample_offspring!(
                forecast,
                day,
                TP_local_sim,
                TP_import_sim,
                TP_indices,
                susceptible_depletion_sim,
                sim,
            )

            # see if the sim is good and break if not
            if day > 0
                reinitialise_allowed =  day >= 7 &&
                    day <= T_observed &&
                    n_restarts < max_restarts &&
                    day - last_injection >= 3

                (bad_sim, injected_cases) = check_sim!(
                    forecast,
                    forecast_start_date,
                    omicron_dominant_date,
                    case_counts,
                    sim,
                    local_cases,
                    window_lengths,
                    min_cases,
                    max_cases,
                    reinitialise_allowed;
                    day=day
                )
            end

            # reset day counter if we injected some cases
            if !injected_cases
                day += 1
            else
                day = start_day[sim]
                n_restarts += 1
                injected_cases = false
                last_injection = day
                if n_restarts > max_restarts
                    bad_sim = true
                end
            end

            # break early if sim inconsistent
            bad_sim && break
        end

        if !bad_sim
            reinitialise_allowed = false
            (bad_sim, injected_cases) = check_sim!(
                forecast,
                forecast_start_date,
                omicron_dominant_date,
                case_counts,
                sim,
                local_cases,
                window_lengths,
                min_cases,
                max_cases,
                reinitialise_allowed,
            )
        end

        if bad_sim
            good_sims[sim] = false
        else
            good_TPs[TP_ind] = true
        end
    end

    # print out the number of good simulations
    println("======================")
    println("Summary for ", state, ":")
    println("- ", sum(good_sims), " good simulations.")
    println("- ", sum(good_TPs), " unique TPs sampled.")
    println("======================")

    # keep only good sim results and truncate the TP to the same period as per D and U
    TP_local_good_truncated = TP_local[(end-T_end):end, good_TPs]
    D_good = forecast.sim_realisations.D[:,:,good_sims]
    U_good = forecast.sim_realisations.U[:,:,good_sims]

    return (D_good, U_good, TP_local_good_truncated)
end


text:
module FEFunctionsTests

using Gridap

using Test
using Gridap

using Gridap.Triangulations: IndexCellFieldWithTriangulation

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))
order = 1
diritag = "boundary"
fespace = H1ConformingFESpace(Float64,model,order,diritag)

trian = Triangulation(model)

ufun(x) = x[1]

uh = interpolate(fespace,ufun)

@test string(uh) == "FEFunction object"
sr = "FEFunction object:\n physdim: 2\n refdim: 2\n valuetype: Float64\n doftype: Float64\n nfree: 9\n ndiri: 16\n ncells: 16"
@test sprint(show,"text/plain",uh) == sr

@test isa(Triangulation(uh),Triangulation)

tags = [7,6]
btrian = BoundaryTriangulation(model,tags)
bquad = CellQuadrature(btrian,degree=0)

buh = restrict(uh,btrian)

cn = integrate(buh,btrian,bquad)
@test isa(cn,CellNumber)
_ = collect(cn)

cf = IndexCellFieldWithTriangulation(uh.cellfield,trian)

quad = CellQuadrature(trian,degree=2)
q = coordinates(quad)
v = collect(evaluate(cf,q))
g = collect(evaluate(gradient(cf),q))

@test isa(∇(cf),IndexCellFieldWithTriangulation)

test_index_cell_field(cf,q,v,g)

@test isa(∇(uh),IndexCellFieldWithTriangulation)
@test isa(Triangulation(∇(uh)),Triangulation)

u = CellField(trian,ufun)

@test isa(u,IndexCellFieldWithTriangulation)

e = u - uh
@test isa(e,IndexCellFieldWithTriangulation)

e = ufun - uh
@test isa(e,IndexCellFieldWithTriangulation)

e = uh - ufun
@test isa(e,IndexCellFieldWithTriangulation)

T = VectorValue{2,Float64}
fespace = H1ConformingFESpace(T,model,order,diritag)
ufun(x) = VectorValue(x[2],x[1])
uh = interpolate(fespace,ufun)

@test isa(div(uh),IndexCellFieldWithTriangulation)
@test isa(curl(uh),IndexCellFieldWithTriangulation)

end # module


text:
export apply_qubit_op!, get_transition_amplitude_with_obs

"""
Divide a qubit operator into the hermite and antihermite parts.
"""
divide_real_imag(op::QubitOperator) =
    (op+hermitian_conjugated(op))/2, (op-hermitian_conjugated(op))/2im

"""
Apply a qubit operator op to |state_ket> and fit the result with
circuit * |state_bra>.
The circuit object will be updated on exit.
The squared norm of op * |state_ket>  will be returned.
state0_bra will not be modified.
"""
function apply_qubit_op!(
    op::QubitOperator,
    state_ket::QuantumState,
    circuit::VariationalQuantumCircuit, state0_bra::QuantumState;
    minimizer=mk_scipy_minimize(),
    verbose=true,
    comm=MPI_COMM_WORLD
    )
    her, antiher = divide_real_imag(op)

    function cost(thetas::Vector{Float64})
        update_circuit_param!(circuit, thetas)
        re_ = get_transition_amplitude_with_obs(circuit, state0_bra, her, state_ket)
        im_ = get_transition_amplitude_with_obs(circuit, state0_bra, antiher, state_ket)
        abs2(1.0 - (re_ + im_ * im))
    end

    thetas_init = get_thetas(circuit)
    if comm !== nothing
        MPI.Bcast!(thetas_init, 0, comm)
    end
    opt_thetas = minimizer(cost, thetas_init)
    norm_right = sqrt(get_expectation_value(hermitian_conjugated(op) * op, state_ket))
    if verbose
       println("norm_right",norm_right)
    end

    update_circuit_param!(circuit, opt_thetas)
    re__ = get_transition_amplitude_with_obs(circuit, state0_bra, her, state_ket)
    im__ = get_transition_amplitude_with_obs(circuit, state0_bra, antiher, state_ket)
    z = re__ + im__ * im
    if verbose
        println("Match in apply_qubit_op!: ", z/norm_right)
    end
    return z
end

"""
Compute <state_bra| circuit^+ obs |state_ket>, where obs is a hermite observable.
"""
function get_transition_amplitude_with_obs(
    circuit::VariationalQuantumCircuit, state0_bra::QuantumState,
    op::QubitOperator,
    state_ket::QuantumState)
    state_bra = copy(state0_bra)
    update_quantum_state!(circuit, state_bra)
    return get_transition_amplitude(op, state_bra, state_ket)
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
struct LogLoss{M<:AbstractMatrix,V<:AbstractVector} <: AbstractLoss
  A::M
  b::V

  function (::Type{LogLoss})(A::AbstractMatrix, b::AbstractVector)
    @assert size(A, 1) == length(b) "`A` and `b` must have compatible dimensions"
    T = float(promote_type(eltype(A), eltype(b)))
    A_ = convert(AbstractMatrix{T}, A)
    b_ = convert(AbstractVector{T}, b)

    new{typeof(A_),typeof(b_)}(A_, b_)
  end
end

nsamples(loss::LogLoss)   = size(loss.A, 1)
nfeatures(loss::LogLoss)  = size(loss.A, 2)

function loss!(loss::LogLoss, x::AbstractVector{T}, g::AbstractVector{T}) where T
  val = -loss.b .* (loss.A * x)
  fval = zero(T)

  for idx in 1:length(val)
    temp = exp(val[idx])
    @inbounds fval += isinf(temp) ? val[idx] : log1p(temp)
    @inbounds val[idx] = isinf(temp) ? -loss.b[idx] : -loss.b[idx]*temp/(one(T) + temp)
  end
  mul!(g, transpose(loss.A), val)
  return fval
end