Packagecompile (#53) close #52

* package compiler seems to have worked

* fixing inclusion of transitive dependencies

* small addition

* adapting package compiler to testdata artifact

and adapting a bit more here and there
- rewritten all tests to use @testitem
- reverting to using [extras]/[targets] instead of test/Project.toml because this works better with @testitem in that dependencies from main do not need to be duplicated
- rerun just precompilestatements

* improving justfile

Author: schlichtanders

.gitignore

@@ -5,3 +5,5 @@
 /docs/build/
 Manifest.toml
 tmp/
+/packagecompiler/app
+/packagecompiler/app.backup

Project.toml

@@ -6,7 +6,6 @@ version = "1.0.0-DEV"
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
@@ -26,12 +25,15 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 Missings = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
 NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 RecurrenceAnalysis = "639c3291-70d9-5ea2-8c5b-839eba1ee399"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 TimeseriesSurrogates = "c804724b-8c18-5caa-8579-6025a0767c70"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 YAXArrayBase = "90b8fcef-0c2d-428d-9c56-5f86629e9d14"
@@ -41,7 +43,6 @@ Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99"
 [compat]
 ArchGDAL = "0.10"
 ArgParse = "1"
-CairoMakie = "0.12"
 ConstructionBase = "1"
 DimensionalData = "0.29"
 DiskArrayTools = "0.1"
@@ -54,18 +55,30 @@ GDAL = "1"
 GeoFormatTypes = "0.4"
 Glob = "1"
 KML = "0.2"
-LinearAlgebra = "1.10.0, 1.11.0"
+LinearAlgebra = "1.10"
 LoggingExtras = "1"
 Missings = "1"
 NetCDF = "0.12"
 Proj = "1"
 Rasters = "0.12,0.13"
 RecurrenceAnalysis = "2"
-Statistics = "1"
+Statistics = "1.10"
+TestItems = "1.0.0"
 TimeseriesSurrogates = "2"
 UnicodePlots = "3"
 YAXArrayBase = "0.6, 0.7"
 YAXArrays = "0.5"
 Zarr = "0.9"
 julia = "1.10"
 
+[extras]
+AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
+
+
+[targets]
+test = ["Test", "TestItemRunner", "Pkg", "Random", "AllocCheck", "BenchmarkTools"]

justfile

@@ -0,0 +1,41 @@
+default:
+  just --list
+
+# create precompile statements used for `just packagecompile`
+precompilestatements:
+    julia --project -e 'import Pkg; Pkg.test(julia_args=`--trace-compile=packagecompiler/precompile_statements.jl`)'
+
+# build an app at `packagecompiler/app/RQADeforestation` with PackageCompiler.jl
+packagecompile:
+    #!/usr/bin/env -S julia --project=packagecompiler
+    using PackageCompiler
+    if isdir("packagecompiler/app") 
+        rm("packagecompiler/app.backup", recursive=true, force=true)
+        mv("packagecompiler/app", "packagecompiler/app.backup")
+    end
+    PackageCompiler.create_app(".", "packagecompiler/app"; 
+        precompile_statements_file="packagecompiler/precompile_statements.jl",
+        # see https://github.com/JuliaLang/PackageCompiler.jl/issues/994
+        include_lazy_artifacts=true,
+        # or try the following, which may work because IntelOpenMP and MKL may actually not be needed
+        # include_transitive_dependencies=false,
+    )
+    # lets have an easy check whether this actually worked
+    touch("packagecompiler/app/done")
+
+# test the app `packagecompiler/app/RQADeforestation` with testdata, writing data to `test/tmp/apptestdata`
+testapp:
+    #!/usr/bin/env bash
+    cd test
+    indir="tmp/apptestdata/in"
+    outdir="tmp/apptestdata/out.zarr"
+    julia --project -e '
+        import Pkg: Artifacts.@artifact_str, ensure_artifact_installed
+        ensure_artifact_installed("rqatestdata", "Artifacts.toml")
+        testdatapath = joinpath(artifact"rqatestdata", "RQADeforestationTestData-1.0")
+        testdir = dirname(ARGS[1])
+        rm(testdir, recursive=true, force=true)
+        mkpath(testdir)
+        cp(testdatapath, ARGS[1])
+    ' -- "$indir"
+    ../packagecompiler/app/bin/RQADeforestation --tile E051N018T3 --continent EU --in-dir "$indir" --out-dir "$outdir"

packagecompiler/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"

packagecompiler/precompile_statements.jl

@@ -7,70 +7,15 @@ using YAXArrays
 using Zarr
 using Distributed: myid
 using NetCDF
+using TestItems
 
 export gdalcube, rqatrend
 
+include("metrics.jl")
 include("auxil.jl")
 include("rqatrend.jl")
 include("analysis.jl")  # TODO what is still needed from analysis now that rqatrend is in its own file?
 include("timestats.jl")
-
-
-function main(tiles=["E048N018T3"]; pol="VH", orbit="*", thresh=3.0)
-    indir = "/eodc/products/eodc.eu/S1_CSAR_IWGRDH/SIG0/"
-    continent = "EU"
-    folders = ["V01R01", "V0M2R4", "V1M0R1", "V1M1R1", "V1M1R2"]
-    corruptedfiles = "corrupted_tiles.txt"
-    # TODO save the corrupt files to a txt for investigation
-    for tilefolder in tiles
-
-        filenamelist = [glob("$(sub)/*$(continent)*20M/$(tilefolder)/*$(pol)_$(orbit)*.tif", indir) for sub in folders]
-
-        allfilenames = collect(Iterators.flatten(filenamelist))
-
-
-        relorbits = unique([split(basename(x), "_")[5][2:end] for x in allfilenames])
-        @show relorbits
-        for relorbit in relorbits
-            for y in [2018, 2019, 2020, 2021, 2022, 2023]
-
-                filenames = allfilenames[findall(contains("$(relorbit)_E"), allfilenames)]
-                @time cube = gdalcube(filenames)
-
-                path = joinpath(YAXDefaults.workdir[], "$(tilefolder)_rqatrend_$(pol)_$(relorbit)_thresh_$(thresh)_year_$(y)")
-                @show path
-                ispath(path * ".done") && continue
-                ispath(path * "_zerotimesteps.done") && continue
-
-                tcube = cube[Time=Date(y - 1, 7, 1) .. Date(y + 1, 7, 1)]
-                @show size(cube)
-                @show size(tcube)
-                if size(tcube, Ti) == 0
-                    touch(path * "_zerotimesteps.done")
-                    continue
-                end
-                try
-                    @time rqatrend(tcube; thresh, outpath=path * ".zarr", overwrite=true)
-                catch e
-
-                    if e.captured.ex isa ArchGDAL.GDAL.GDALError
-                        println("Found GDALError:")
-                        println(e.captured.ex.msg)
-                        continue
-                    else
-                        rethrow(e)
-                    end
-                end
-                #=@everywhere begin
-                    fname = "$(VERSION)_$(getpid())_$(time_ns()).heapsnapshot"
-                    Profile.take_heap_snapshot(fname;streaming=true)
-                end
-                =#
-
-                touch(path * ".done")
-            end
-        end
-    end
-end
+include("main.jl")
 
 end

src/RQADeforestation.jl

@@ -120,7 +120,9 @@ function gdalcube(filenames::AbstractVector{<:AbstractString})
     stacked_gdbs = diskstack(all_gdbs)
     attrs = copy(gdb.attrs)
     #attrs["add_offset"] = Float16(attrs["add_offset"])
-    attrs["scale_factor"] = Float16(attrs["scale_factor"])
+    if haskey(attrs, "scale_factor")
+        attrs["scale_factor"] = Float16(attrs["scale_factor"])
+    end
     all_cfs = CFDiskArray(stacked_gdbs, attrs)
     return YAXArray((onecube.axes..., taxis), all_cfs, onecube.properties)
     #datasetgroups = [datasets[group] for group in groupinds]

src/auxil.jl

@@ -0,0 +1,130 @@
+using ArgParse
+using YAXArrays: YAXDefaults
+
+const argparsesettings = ArgParseSettings()
+@add_arg_table! argparsesettings begin
+    "--threshold", "-t"
+    help = "Threshold for the recurrence matrix computation"
+    default = 3.0
+
+    "--polarisation", "-p"
+    help = "Polarisation that should be stacked"
+    default = "VH"
+
+    "--years", "--year", "-y"
+    help = "Year in which the RQA Trend should be detected. 
+    We take a buffer of six month before and after the year to end up with two years of data."
+    default = [2018, 2019, 2020, 2021, 2022, 2023]
+    nargs = '+'
+    arg_type = Int
+
+    "--orbit", "-o"
+    help = "One of: Orbit number, 'A' for ascending, 'D' for descending, '*' for all orbits"
+    default = "*"
+
+    "--out-dir", "-d"
+    help = "Path to output zarr dataset"
+    default = "out.zarr"
+    dest_name = "outdir"
+
+    "--in-dir"
+    help = "Path to input"
+    required = true
+    dest_name = "indir"
+
+    "--continent"
+    help = "continent code for the tile to be processed"
+    required = true
+
+    "--tiles", "--tile"
+    help = "Tile id to be processed"
+    required = true
+    nargs = '+'
+    arg_type = String
+
+    "--folders", "--folder"
+    help = "subfolders taken into account"
+    default = ["V01R01", "V0M2R4", "V1M0R1", "V1M1R1", "V1M1R2"]
+    arg_type = String
+    nargs = '*'
+end
+
+
+function julia_main()::Cint
+    main(; Dict(Symbol(k) => v for (k, v) in parse_args(argparsesettings))...)
+    return 0
+end
+
+
+function main(;
+    tiles::Vector{String},
+    continent::String, 
+    indir::String,
+    outdir="out.zarr",
+    years=[2018, 2019, 2020, 2021, 2022, 2023],
+    polarisation="VH", 
+    orbit="*", 
+    threshold=3.0, 
+    folders=["V01R01", "V0M2R4", "V1M0R1", "V1M1R1", "V1M1R2"]
+)
+    if isdir(indir) && isempty(indir)
+        error("Input directory $indir must not be empty")
+    end
+    if isdir(outdir)
+        @warn "Resume from existing output directory"
+    else
+        mkdir(outdir)
+        @info "Write output to $outdir"
+    end
+    YAXDefaults.workdir[] = outdir
+
+    corruptedfiles = "corrupted_tiles.txt"
+    # TODO save the corrupt files to a txt for investigation
+    for tilefolder in tiles
+
+        filenamelist = [glob("$(sub)/*$(continent)*20M/$(tilefolder)/*$(polarisation)_$(orbit)*.tif", indir) for sub in folders]
+        allfilenames = collect(Iterators.flatten(filenamelist))
+        @show allfilenames
+
+        relorbits = unique([split(basename(x), "_")[5][2:end] for x in allfilenames])
+        @show relorbits
+        for relorbit in relorbits
+            for y in years
+
+                filenames = allfilenames[findall(contains("$(relorbit)_E"), allfilenames)]
+                @time cube = gdalcube(filenames)
+
+                path = joinpath(YAXDefaults.workdir[], "$(tilefolder)_rqatrend_$(polarisation)_$(relorbit)_thresh_$(threshold)_year_$(y)")
+                @show path
+                ispath(path * ".done") && continue
+                ispath(path * "_zerotimesteps.done") && continue
+
+                tcube = cube[Time=Date(y - 1, 7, 1) .. Date(y + 1, 7, 1)]
+                @show size(cube)
+                @show size(tcube)
+                if size(tcube, 3) == 0
+                    touch(path * "_zerotimesteps.done")
+                    continue
+                end
+                try
+                    @time rqatrend(tcube; thresh=threshold, outpath=path * ".zarr", overwrite=true)
+                catch e
+
+                    if hasproperty(e, :captured) && e.captured.ex isa ArchGDAL.GDAL.GDALError
+                        println("Found GDALError:")
+                        println(e.captured.ex.msg)
+                        continue
+                    else
+                        rethrow(e)
+                    end
+                end
+                #=@everywhere begin
+                    fname = "$(VERSION)_$(getpid())_$(time_ns()).heapsnapshot"
+                    Profile.take_heap_snapshot(fname;streaming=true)
+                end
+                =#
+                touch(path * ".done")
+            end
+        end
+    end
+end

src/main.jl

@@ -0,0 +1,14 @@
+import Distances
+struct CheckedEuclidean <: Distances.UnionMinkowskiMetric end
+
+@inline (dist::CheckedEuclidean)(a, b) = Distances._evaluate(dist, a, b, Distances.parameters(dist))
+@inline Distances.eval_op(::CheckedEuclidean, ai::Number, bi::Number) = abs(ai - bi)
+@inline Distances.eval_op(::CheckedEuclidean, ai::Integer, bi::Integer) = Base.Checked.abs(Base.Checked.checked_sub(ai, bi))
+Distances.eval_end(::CheckedEuclidean, s) = s
+
+
+@testitem "CheckedEuclidean" begin
+    # this nested using seems necesarry because of using test/Project.toml
+    using Distances
+    @test evaluate(RQADeforestation.CheckedEuclidean(), 2, 5) == evaluate(Euclidean(), 2, 5)
+end