{"cells":[{"metadata":{},"cell_type":"markdown","source":"# COLLABORATION WELCOME!\n***If the ideas, code, features, or predictions below are interesting or helpful to you, please reach out! Given how different this appraoch is to the common Python + Deep Learning approach, there might be opportunity for us to team up and ensemble our models. I'm releasing my code and ideas publicly here, but if you use significant parts of this to make your solution better, please consider inviting me to your team ^_^***"},{"metadata":{},"cell_type":"markdown","source":"## CHANGELOG\n\n### Update 2021.03.27\n\nI've overhauled the notebook to be aimed more at an actual submission.\n- It can run offline, thanks to the downloads being moved to another notebook.\n- I've implemented some hefty feature engineering, including spatially-aware features, in Julia.\n- I've moved the modeling over to LightGBM in Python.\n- I show off my RLE encoding/decoding skills to move segmentation masks from Julia to Python (I use numba to JIT compile Python code to get a 100x speedup for reading RLE masks)"},{"metadata":{},"cell_type":"markdown","source":"## Introduction\n\nThis notebook is a bit atypical for a Kaggle image competition. First, it's written purely in Julia, not Python. Second, it uses no deep learning. If this sounds interesting to you, please continue reading!\n\n### Goals\n\nThe primary goal of this notebook is to take a crack at the competition from scratch (segmentation and prediction) with a an atypical set of methods and a different programming language and package ecosystem. The secondary goal is to try and run everything fast on a CPU. \n\n### Backstory\n\nI'm new to the Julia programming language, but I'm very excited about it. I got interested in this competition from a friend just after wrapping up my first project in Julia, and I thought to myself, \"hey, let's take a crack at doing this in Julia, just as a learning exercise.\" Since it's clunky actually developing Julia in Kaggle notebooks (no autocomplete, can't interrup long-running computation, etc.), I started working locally on my laptop, which means no GPU. It's harder (but possible!) to write Julia code for GPU than CPU, anyhow, and as I'm just starting to figure things out, I figured I'd try to make my approach efficient enough to run on a laptop CPU. Thus the goals above came about naturally."},{"metadata":{},"cell_type":"markdown","source":"## Installing Julia\nBefore we can run any code, we have to install the Julia language. Sadly, even though Julia is the \"Ju\" in **Ju**pyter, Kaggle no longer supports Julia Kernels directly, so we have to use Julia through a Python kernel, which is a bit painful.\n\nIn order to make this notebook work offline, I separated most of the installation steps into [another notebook, which you can take a look at if you're curious](https://www.kaggle.com/lukemerrick/julia-download)."},{"metadata":{"trusted":true},"cell_type":"code","source":"!cd ../input/julia-download/install_julia && sh install.sh","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install --no-index --find-links ../input/pycocotools/ pycocotools","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nfrom multiprocessing import cpu_count\nos.environ[\"JULIA_NUM_THREADS\"] = str(cpu_count())\nimport julia\nfrom julia.api import Julia\n# cannot use precompiled packages with pyjulia on linux :-(\n#     but we can use a system image with PyCall pre-compiled to speed that up\njl = Julia(\n    sysimage=\"../input/julia-download/install_julia/python_julia_sysimage.so\",\n    compiled_modules=False\n)\n%load_ext julia.magic","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%julia\n# confirm we're using more than one thread\nusing Base.Threads: nthreads\nnthreads()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## \"Installing\" home-baked dependencies\nIn addition to making a semi-polished public package, I also created modules to organize utility functions. Unfortunately this buries the logic for segmentation a bit. I'm planning to come back and run the lower-level functions and better explain the segmentation details, but for now, I'm working to getting things going end-to-end."},{"metadata":{"trusted":true},"cell_type":"code","source":"!mkdir ProteinAtlas","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"%%writefile ProteinAtlas/ProteinAtlas.jl\nmodule ProteinAtlas\n\nimport\n    CellSegmentation,\n    FileIO,\n    ImageContrastAdjustment,\n    ImageCore,\n    ImageSegmentation\n\nusing ColorTypes: Gray, HSV, N0f8, RGB\nusing CellSegmentation: CellImage, CellImageNoAxes, cell_image, nearest_neighbor_resize\nusing ImageTransformations: imresize\nusing Statistics: quantile, median, mean\n\nexport\n    HPA_CHANNELS,\n    load_example,\n    rgb_hpa_image,\n    HPASegmentation,\n    segment_cells,\n    SegmentedCell,\n    segmentation_image,\n    extract_cell,\n    resize_cell\n\n\nconst HPA_CHANNELS = [:mt, :er, :nu, :pr]\n\n\"\"\"\nLoads all four channels of an example as Grayscale images.\n* red microtubule \"mt\" channel\n* yellow endoplasmic reticulum \"er\" channel\n* blue nucleus \"nu\" channel\n* green protein of interest \"pr\" channel\n\nThen offers a StackedView of them together.\n\"\"\"\nfunction load_example(dir::String, key::String)::CellImage\n    mt = FileIO.load(joinpath(dir, \"$(key)_red.png\"))\n    er = FileIO.load(joinpath(dir, \"$(key)_yellow.png\"))\n    nu = FileIO.load(joinpath(dir, \"$(key)_blue.png\"))\n    pr = FileIO.load(joinpath(dir, \"$(key)_green.png\"))\n    example_data = Array{Gray{N0f8},3}(\n        undef, (size(mt)..., 4)\n    )\n    example_data[:, :, 1] = mt\n    example_data[:, :, 2] = er\n    example_data[:, :, 3] = nu\n    example_data[:, :, 4] = pr\n    return cell_image(example_data, HPA_CHANNELS)\nend\n\n\"\"\"\nConvert a HPAImage to RGB for visualization.\n    * red for the microtubule channel\n    * yellow for the endoplasmic reticulum channel\n    * blue for the nucleus channel\n    * green for the protein channel\n\nOptions\n-------\n`include_er`: Toggles the yellow er channel, so you can have true RGB of mu, pr, nu\n`match_intensity_histograms`: Scales the brightness of each channel to match the average\n    across channels. Generally makes the image more boring, but also easier to\n    see locations of all three channels (no one channel should dominate).\n`gamma_adjust`: Gamma adjustment. Between dropping to 0.5-0.8 can brighten\n    dark stuff. Especially useful when match_intensity_histograms=true.\n`brightness_scale`: Linearly scale pixel values up by this amount\n    (e.g. 1.5 is a 50% increase in birghtness)\n\"\"\"\nfunction rgb_hpa_image(\n    image::Union{CellImage,CellImageNoAxes};\n    include_er::Bool=true,\n    match_intensity_histograms::Bool=false,\n    gamma_adjust::Real=1.0,\n    brightness_scale::Real=1.0\n)::Array{RGB{N0f8},2}\n    # restore channel axes if missing\n    if isa(image, CellImageNoAxes)\n        image = cell_image(image, HPA_CHANNELS)\n    end\n\n    # drop :er channel if we're skipping it\n    if !include_er\n        image = image[channel=[:mt, :pr, :nu]]\n    end\n\n    # optionally rescale each channel so its color histogram matches that of\n    # the average across all channels\n    if match_intensity_histograms\n        image = deepcopy(image)\n        channel_avg = mean(image; dims=3)[:, :, 1]\n        adjustment = ImageContrastAdjustment.Matching(; targetimg=channel_avg)\n        for channel_slice in eachslice(image; dims=3)\n            ImageContrastAdjustment.adjust_histogram!(channel_slice, adjustment)\n        end\n    end\n\n    # red, green, and blue are channels directly\n    red = image[channel=:mt]\n    green = image[channel=:pr]\n    blue = image[channel=:nu]\n\n    # we treat yellow is as an even mix of red and green by adding to both channels\n    if include_er\n        red += image[channel=:er] / 2\n        green += + image[channel=:er] / 2\n        ImageCore.clamp01!(red)\n        ImageCore.clamp01!(green)\n    end\n    result = RGB{N0f8}.(red, green, blue)\n    gamma_correction = ImageContrastAdjustment.GammaCorrection(gamma=gamma_adjust)\n    ImageContrastAdjustment.adjust_histogram!(result, gamma_correction)\n    result = RGB{N0f8}.(ImageCore.clamp01.(brightness_scale .* result))\n    return result\nend\n\nstruct HPASegmentation\n    nuclei_segmentation_map::Array{Int64,2}\n    cell_segmentation_map::Array{Int64,2}\nend\n\n\"\"\"\nSegments nuclei and uses that to segment the whole cells.\n\nAlgorithm used is marker-based watershed.\n\nOptions\n-------\n`resize_px`: Internally the image is downscaled to resize_px x resize_px for performance.\n    Smaller is faster, but too small may give poor results. \n`nucleus_bg_arguments`: Arguments passed to `get_channel_background` for computing\n    nucleus background.\n`cell_bg_arguments`: Arguments passed to `get_channel_background` for computing\n    whole cell background.\n`nucleus_marker_quantile`: Quantile of nonzero distance values to use in\n    identifying nucleus segmentation markers. Higher will result in\n    smaller/fewer markers. Too high and small/noisy nuclei will be missed, too\n    low and close nuclei may merge.\n\"\"\"\nfunction segment_cells(\n    image::CellImage;\n    resize_px::Int64=512,\n    nucleus_bg_arguments::Dict{Symbol,Float64}=Dict(\n        :median_filter_window_size => 0.005,\n        :background_nonzero_quantile => 0.2,\n        :area_opening_window_size => 0.02\n    ),\n    cell_bg_arguments::Dict{Symbol,Float64}=Dict(\n        :median_filter_window_size => 0.01,\n        :background_nonzero_quantile => 0.2,\n        :area_opening_window_size => 0.02\n    ),\n    nucleus_marker_quantile::Float64=0.4,\n)::HPASegmentation\n    # downscale image \n    small_image = cell_image(imresize(image, (resize_px, resize_px)), HPA_CHANNELS)\n    # segment nuclei using marker-based watershep\n    nucleus_background = CellSegmentation.get_background_mask(\n        small_image[channel=:nu].data; nucleus_bg_arguments...\n    )\n    nucleus_seg = CellSegmentation.segment_via_background(\n        nucleus_background; watershed_marker_quantile=nucleus_marker_quantile\n    )\n    nucleus_seg_map = ImageSegmentation.labels_map(nucleus_seg)\n\n    # segment full cells\n    cell_grayscale = (\n        small_image[channel=:nu].data\n        + small_image[channel=:er].data\n        + small_image[channel=:mt].data\n    )\n    cell_background = CellSegmentation.get_background_mask(\n        cell_grayscale; cell_bg_arguments...\n    )\n    cell_seg = CellSegmentation.segment_via_background(\n        cell_background; markers=nucleus_seg_map\n    )\n    cell_seg_map = ImageSegmentation.labels_map(cell_seg)\n\n    # upscale segmentation maps\n    full_size = size(image)[1:2]\n    nucleus_seg_map = Int.(nearest_neighbor_resize(nucleus_seg_map, full_size))\n    cell_seg_map = Int.(nearest_neighbor_resize(cell_seg_map, full_size))\n    hpa_segmentation = HPASegmentation(nucleus_seg_map, cell_seg_map)\nend\n\n# create evenly spaced colors palette using HSV across various brightness levels\nn_angles = 8\nn_saturations = 3\ncolors = [\n    HSV(θ + (i % 2) * 360/n_angles/2, 0.5 + i / n_saturations / 2, 0.7)\n    for θ in 360/n_angles:360/n_angles:360, i in n_saturations:-1:1\n][:]\npushfirst!(colors, HSV(0, 0, 0))\n\n# define function to max out saturation of nuclei pixels\nmax_brightness(color::HSV)::HSV = HSV(color.h, color.v, 1)\n\n\"\"\"\nCreate HSV image of a segmentation map\n\"\"\"\nfunction segmentation_image(seg::HPASegmentation)::AbstractMatrix{HSV}\n    seg_image = map(i -> colors[i % length(colors) + 1], seg.cell_segmentation_map)\n    nucleus_mask = seg.nuclei_segmentation_map .> 0\n    seg_image[nucleus_mask] = max_brightness.(seg_image[nucleus_mask])\n    return seg_image\nend\n\nstruct SegmentedCell\n    image::CellImage\n    nucleus_mask::BitMatrix\n    cell_mask::BitMatrix\nend\n\n\"\"\"\nExtract the contents of a segment mask and place it on a black background.\n\nAlso returns indices of pixesl in the nucleus and indices of pixels in the whole cell.\n\"\"\"\nfunction extract_cell(\n    nucleus_mask::AbstractArray,\n    cell_mask::AbstractArray,\n    full_image::CellImage\n)::SegmentedCell\n    nucleus_indices = findall(nucleus_mask)\n    cell_indices = findall(cell_mask)\n    idx_min = minimum(cell_indices)\n    idx_max = maximum(cell_indices)\n    shift = idx_min - CartesianIndex(1, 1)\n    shifted_cell_indices = cell_indices .- (shift,)\n    shifted_nucleus_indices = nucleus_indices .- (shift,)\n    max_dim = maximum((idx_max - idx_min).I) + 1\n    channel_names = full_image.axes[3].val\n    extracted_image = zeros(Gray{N0f8}, (max_dim, max_dim, length(channel_names)))\n    extracted_image = cell_image(extracted_image, channel_names)\n    extracted_image[shifted_cell_indices, :] = full_image[cell_indices, :]\n    out_nucleus_mask = falses(max_dim, max_dim)\n    out_cell_mask = falses(max_dim, max_dim)\n    out_nucleus_mask[shifted_nucleus_indices] .= 1\n    out_cell_mask[shifted_cell_indices] .= 1\n    return SegmentedCell(extracted_image, out_nucleus_mask, out_cell_mask)\nend\n\nfunction resize_cell(cell::SegmentedCell, new_size)\n    return SegmentedCell(\n        CellSegmentation.cell_image(imresize(cell.image, new_size), HPA_CHANNELS),\n        Bool.(nearest_neighbor_resize(cell.nucleus_mask, new_size)),\n        Bool.(nearest_neighbor_resize(cell.cell_mask, new_size))\n    )\nend\n\nend","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# TODO: Walk through individual components\nIn this latest rounds of updates, I'm pushing to getting my work running end-to-end (including submission of predictions). In future updates, I plan to come back and show step-by-step how the segmentation and feature engineering work."},{"metadata":{},"cell_type":"markdown","source":"# Putting it all together (at least the Julia part)\nBelow is a script that handles the loading, segmentation, and feature-encoding of all the cell images. As output, we get feature vectors for all the cells identified by segmenting all the train and test images, as well as JSON files which contain RLE-encoded segmentation masks for all the test cells."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n%%julia\n@info \"Running Julia imports. Compiling can take a few minutes, so be patient.\"\nimport\n    CellSegmentation,\n    CSV,\n    FileIO,\n    ImageMorphology,\n    ImageSegmentation,\n    JSON,\n    Parquet,\n    Random\n\nusing Base.Threads: @threads\nusing CellSegmentation: nearest_neighbor_resize, CellImage\nusing ColorTypes: RGB, N0f8, HSV, Gray\nusing DataFrames\nusing ImageCore: clamp01\nusing ImageTransformations: imresize\nusing MosaicViews: mosaicview\nusing ProgressMeter: Progress, next!, @showprogress\nusing Statistics: quantile, std, mean, cor\nusing StatsBase: counts, rle\n\n# local modules\npush!(LOAD_PATH, \"./ProteinAtlas\")\nusing ProteinAtlas","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n%%julia\nconst N_TRAINING_IMG = 20\n@info \"Starting main process! JIT compiling can take a long time, so be patient.\"\n@info \"Subsetting to $N_TRAINING_IMG images.\"\n\nconst TRAIN_DIR = \"../input/hpa-single-cell-image-classification/train\"\nconst TEST_DIR = \"../input/hpa-single-cell-image-classification/train\"\nconst TRAIN_CSV = \"../input/hpa-single-cell-image-classification/train.csv\"\n\nstruct RLEMask\n    mask_nrow::Integer\n    mask_ncol::Integer\n    run_values::Vector{<:Real}\n    run_lengths::Vector{<:Integer}\n    function RLEMask(mask::BitMatrix)\n        mask_nrow, mask_ncol = size(mask)\n        run_values, run_lengths = rle(mask[:])\n        return new(mask_nrow, mask_ncol, UInt8.(run_values), run_lengths)\n    end\nend\n\n\nfunction sample_rows(df; n=10, seed=0)\n    rng = Random.MersenneTwister(seed)\n    n = min(n, nrow(df))\n    indices = Random.randperm(nrow(df))[1:n]\n    return df[indices, :]\nend\n\n\"\"\"\nSplit up a blob mask (nucleus or whole cell) into rings and pie slices\n\"\"\"\nfunction region_segment(mask:: AbstractMatrix; n_ring = 4, n_angle=8)\n    height, width = size(mask)\n    bg_mask = .!mask\n    distance_from_bg = ImageMorphology.distance_transform(\n        ImageMorphology.feature_transform(bg_mask)\n    )\n    lo, hi = extrema(distance_from_bg[:])\n    slice_width = (hi - lo) / n_ring\n    ring_masks = falses(height, width, n_ring)\n    for i in 1:n_ring\n        ring_masks[:, :, i] = (\n            lo + slice_width * (i - 1) .< distance_from_bg .<= lo + slice_width * i\n        )\n    end\n    \n    # angle slices\n    _, center = findmax(distance_from_bg)\n    center_r, center_c = center.I\n    offset_min = CartesianIndex(1 - center_r, 1 - center_c)\n    offset_max = CartesianIndex(height - center_r, width - center_c)\n    polar_angles = [atan(x.I...) for x in offset_min:offset_max]\n    slice_angle = 2 * π / n_angle\n    slice_masks = falses(height, width, n_angle)\n    for i in 1:n_angle\n        slice = @. ((i - 1) * slice_angle - π <= polar_angles < i * slice_angle - π)\n        slice_masks[:, :, i] = slice .& mask\n    end\n    return ring_masks, slice_masks\nend\n\n\n\"\"\"\nSegment the bright \"speckles\" in the protein channel\n\"\"\"\nfunction segment_bright_protein_shapes(protein_channel; bright_quantile=0.98)\n    nonzero_pr = protein_channel .> 0\n    threshold = any(nonzero_pr) ? quantile(protein_channel[nonzero_pr], bright_quantile) : 0.0\n    protein_mask = protein_channel .> threshold\n    seg = ImageMorphology.label_components(protein_mask)\n    return seg\nend\n\n\n\"\"\"\nExtract a fixed-length vector of numerical summary statistics for a region\n    - min, max, median, iqr, mean, std of each color channel\n    - bi-channel correlation for each color channel\n    - protein speckles per pixel\n    - protein speckle size distribution: max, median, iqr, mean, std\n\"\"\"\nfunction summarize_region(\n    image::CellImage,\n    region_mask::AbstractMatrix,\n    protein_speckles::AbstractMatrix;\n    aggregations=(minimum, maximum, mean, std),\n    quantiles=(0.05, 0.25, 0.5, 0.75, 0.95)\n)::Vector{Float32}\n    result = Float32[]\n    region_pixels = Float32.(image[region_mask, :])\n    region_pixels = reshape(region_pixels, size(region_pixels)[1:2])  # squeeze\n    sort!(region_pixels; dims=1)\n    segment_speckles = protein_speckles[region_mask]\n\n    # univariate summaries of pixel information\n    for agg_fn in aggregations\n        append!(result, agg_fn(region_pixels; dims=1)[:])\n    end\n\n    # quantiles\n    for channel_pixels in eachslice(region_pixels; dims=2)\n        append!(result, quantile(channel_pixels, quantiles; sorted=true))\n    end\n\n    # bi-channel correlations\n    pixel_correlations = cor(region_pixels; dims=1)[\n        [CartesianIndex(i, j)\n        for i in Base.OneTo(size(region_pixels, 2))\n            for j in Base.OneTo(i - 1)]\n    ]\n    append!(result, pixel_correlations)\n\n    # protein speckle information\n    speckle_sizes = counts(segment_speckles)[2:end]  # drop background (which is value 0)\n    speckle_sizes = speckle_sizes[speckle_sizes .!= 0]  # drop empty\n    if length(speckle_sizes) == 0\n        speckle_frac = 0\n        speckle_frac_per_pixel = 0\n        speckle_aggregations = zeros(length(aggregations))\n        speckle_quantiles = zeros(length(quantiles))\n        speckles_per_pixel = 0\n    else\n        n_total_speckles = maximum(protein_speckles)\n        n_pixels_in_region = length(segment_speckles)\n        # fraction of total speckles in this region\n        speckle_frac = length(speckle_sizes) / n_total_speckles\n        speckle_frac_per_pixel = speckle_frac / n_pixels_in_region\n        speckle_aggregations = [f(speckle_sizes) for f in aggregations]\n        speckle_quantiles = quantile(speckle_sizes, quantiles)\n        speckles_per_pixel = length(speckle_sizes) / n_pixels_in_region\n    end\n    push!(result, speckle_frac)\n    push!(result, speckle_frac_per_pixel)\n    append!(result, speckle_aggregations)\n    append!(result, speckle_quantiles)\n    push!(result, speckles_per_pixel)\n    replace(result, NaN => 0)\n    return result\nend\n\nfill_na_div(a, b) = replace!(a ./ b, NaN => 1, Inf => 10, -Inf => -10)\n\n\n\"\"\"\nExtracts a fixed-length vector of numerical summary statistics for a cell_image\n    - region-level statistics for cell, nucleus, and nonnucleus\n    - (max - min) range of region-level statistics across ring and slice regions of the nucleus\n    - (max - min) range of region-level statistics across ring and slice regions of the cell\n    - ratios of region-level statistics and ranges \n\"\"\"\nfunction summarize_cell(cell::SegmentedCell)\n    features = Float32[]\n    protein_speckles = segment_bright_protein_shapes(cell.image[channel=:pr])\n    cell_summary = summarize_region(cell.image, cell.cell_mask, protein_speckles)\n    nucleus_summary = summarize_region(cell.image, cell.nucleus_mask, protein_speckles)\n    nonnucleus_mask = cell.cell_mask .& .!cell.nucleus_mask\n    nonnucleus_summary = (\n        any(nonnucleus_mask) ?\n        summarize_region(cell.image, nonnucleus_mask, protein_speckles) :\n        cell_summary  # copy the cell/nucleus summary if no nonnucleus region exists\n    )\n    nucleus_vs_nonnucleus = fill_na_div(nucleus_summary, nonnucleus_summary)\n    append!(features, cell_summary)\n    append!(features, nucleus_summary)\n    append!(features, nonnucleus_summary)\n    append!(features, nucleus_vs_nonnucleus)\n    for ((ring_masks, slice_masks), summary) in (\n        (region_segment(cell.cell_mask), cell_summary),\n        (region_segment(cell.nucleus_mask), nucleus_summary),\n    )\n        ring_summaries = cat(\n            (\n                summarize_region(cell.image, region_mask, protein_speckles)\n                for region_mask in eachslice(ring_masks, dims=3)\n                if sum(region_mask) > 0\n            )...;\n            dims=2\n        )\n        slice_summaries = cat(\n            (\n                summarize_region(cell.image, region_mask, protein_speckles)\n                for region_mask in eachslice(slice_masks, dims=3)\n                if sum(region_mask) > 0  # happens when cell on edge of image\n            )...;\n            dims=2\n        )\n        ring_range = [max - min for (min, max) in extrema(ring_summaries; dims=2)]\n        slice_range = [max - min for (min, max) in extrema(slice_summaries; dims=2)]\n        normalized_ring_range = fill_na_div(ring_range, summary)\n        normalized_slice_range = fill_na_div(slice_range, summary)\n        ring_vs_slice = fill_na_div(ring_range, slice_range)\n        append!(features, normalized_ring_range)\n        append!(features, normalized_slice_range)\n        append!(features, ring_vs_slice)\n    end\n    return features\nend\n\nfunction load_and_featurize(dir, key)\n    feature_vectors = Vector{Float32}[]\n    rle_masks = RLEMask[]\n    image = load_example(dir, key)\n    seg = segment_cells(image; resize_px=400)\n    n_cells = maximum(seg.cell_segmentation_map)\n    for cell_number in 1:(n_cells - 1)\n        nucleus_mask = seg.nuclei_segmentation_map .== cell_number\n        cell_mask = seg.cell_segmentation_map .== cell_number\n        cell = extract_cell(nucleus_mask, cell_mask, image)\n        cell = resize_cell(cell, min(size(cell.nucleus_mask), resize_size))\n        try\n            push!(feature_vectors, summarize_cell(cell))\n            push!(rle_masks, RLEMask(cell_mask))  # only push the RLE mask if featurization succeeds\n        catch\n            @error \"Error featurizing cell $cell_number of image $key\"\n        end\n    end\n    feature_matrix = nothing\n    if length(feature_vectors) > 0\n        feature_matrix = mapreduce(transpose, vcat, feature_vectors)\n    else\n        @error \"No feature vectors for image $key\"\n    end\n    return feature_matrix, rle_masks\nend\n\n\n#####\n##### Key-label data loading and sampling\n#####\n\n# load train.csv and parse/multihot-encode labels resulting in columns :id, :0, :1, ...\n@info \"Loading train.csv\"\ntrain = CSV.read(TRAIN_CSV, DataFrame)\nparse_labels(column) = [parse(Int, x) for x in split(column, \"|\")]\nselect!(train, :ID => :id, :Label => ByRow(parse_labels) => :label)\nsample_train = (\n    N_TRAINING_IMG === nothing ?\n    train :\n    sample_rows(train, n=N_TRAINING_IMG)\n\n)\n\n\n#####\n##### Training image loading/segmentation/featurization\n#####\n@info \"[Train] Loading, segmenting, and featurizing cells from $(nrow(sample_train)) images\"\nresize_size = (200, 200)\nprogress = Progress(nrow(sample_train); desc=\"Images...\")\nres = Vector(undef, nrow(sample_train))\nenumerated_itr = collect(enumerate(sample_train[!, :id]))\n\n# precompile trigger\nload_and_featurize(TRAIN_DIR, enumerated_itr[1][2])\n\n# parallel processing of load/segment/extract/featurize\n@threads for (i, key) in enumerated_itr\n    feature_matrix, rle_masks = load_and_featurize(TRAIN_DIR, key)\n    res[i] = (key, feature_matrix)\n    next!(progress)\nend\n\n# remove empty results\n@info \"Removing empty results for images with no cells found\"\ni = 1\nwhile i <= length(res)\n    if res[i][2] === nothing\n        popat!(res, i)\n        global i = 1\n    end\n    i += 1\nend\n\n# create dataframe of features giant list of masks\n@info \"Concatenating results into feature table\"\nfeature_df_list = DataFrame[]\nfor (key, feature_matrix) in res\n    df = DataFrame(feature_matrix)\n    insertcols!(df, 1, :key => key)\n    push!(feature_df_list, df)\nend\nfeature_df = vcat(feature_df_list...)\n\n@info \"Writing features to parquet\"\nParquet.write_parquet(\"train_features.parquet\", feature_df)\n\n\n#####\n##### Test image loading/segmentation/featurization\n#####\ntest_image_keys = unique([\n    String(split(x, \"_\")[1])\n    for x in readdir(TEST_DIR)\n    if endswith(x, \".png\")\n])\nif N_TRAINING_IMG !== nothing\n    test_image_keys = test_image_keys[1:N_TRAINING_IMG]\nend\n\n@info \"[Test] Loading, segmenting, and featurizing cells from $(length(test_image_keys)) images\"\nprogress = Progress(length(test_image_keys); desc=\"Images...\")\nres = Vector(undef, length(test_image_keys))\nenumerated_itr = collect(enumerate(test_image_keys))\n\n# parallel processing of load/segment/extract/featurize\n@threads for (i, key) in enumerated_itr\n    feature_matrix, rle_masks = load_and_featurize(TEST_DIR, key)\n    res[i] = (key, feature_matrix, rle_masks)\n    next!(progress)\nend\n\n# remove empty results\n@info \"Removing empty results for images with no cells found\"\ni = 1\nwhile i <= length(res)\n    if res[i][2] === nothing\n        popat!(res, i)\n        global i = 1\n    end\n    i += 1\nend\n\n# create dataframe of features giant list of masks\n@info \"Concatenating results into feature table and RLE-encoded mask list\"\nfeature_df_list = DataFrame[]\nrle_mask_list = RLEMask[]\nfor (key, feature_matrix, rle_masks) in res\n    df = DataFrame(feature_matrix)\n    insertcols!(df, 1, :key => key)\n    push!(feature_df_list, df)\n    append!(rle_mask_list, rle_masks)\nend\nfeature_df = vcat(feature_df_list...)\n\n@info \"Writing features to parquet\"\nParquet.write_parquet(\"test_features.parquet\", feature_df)\n\n@info \"Writing rle masks to JSON\"\n# JSON requires the entire file to be read at once, so we need to chunk it\n#  if we want the reader to be able to read some at a time\nout_dir = \"rle_masks\"\nrm(out_dir; force=true)\nmkpath(out_dir)\nmasks_per_file = 5_000\nfor (file_i, mask_i) in enumerate(1:masks_per_file:length(rle_mask_list))\n    end_i = min(mask_i + masks_per_file - 1, length(rle_mask_list))\n    write(joinpath(out_dir, \"$file_i.json\"), JSON.json(rle_mask_list[mask_i:end_i]))\nend\n\nprintln(\"Done!\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Modeling in Python\n\nI've decided to switch over to Python for the modeling section here at the end.\n\nWhy? I have now completed Julia code to segment the images and compute spatially-aware statistics that encode each cell as a fixed-length numeric feature vector, and the last step of fitting a model is relatively small, after all. However, I have started losing steam on this project, and doing the remainder in Julia would take much more work than just using Python.\n\nFor one, I'm familiar with several popular ML libraries in Python, but I'm not yet familiar with any of the Julia ML ecosystem. In a previous version of this notebook, I just used a Newton-Raphson Logistic Regression solver I wrote from scratch in Julia. Unfortunately, with my more serious feature enginering code I now have over 500 features per image, about 30x as many features as before. The Hessian inversion in the full Newton-Raphson solver could cost $30^3 = 27,000$ times the compute, and so it doesn't seem feasible to continue using that solver.\n\nAdditionally, it seems pretty nonsensical to try and re-implement in Julia the complicated segmentation mask encoding function needed for submission. Instead, it's much easier to export a more standard RLE encoding of the masks into JSON from Julia, and then to use Python to re-encode via the function provided by the organizers and the `pycocotools` library.\n\n\nThe Python code below basically goes through the typical GBM baseline approach to non-computer-vision Kaggle competitions. We load the features into DataFrames, feed that into LightGBM, and boom, we get predictions. There is a little nuance in that we train one model per class independently (a naive approach ignoring class interaction)"},{"metadata":{"trusted":true},"cell_type":"code","source":"import base64\nimport itertools\nimport json\nimport multiprocessing\nimport zlib\nfrom functools import reduce\nfrom pathlib import Path\n\nimport lightgbm as lgb\nimport numpy as np\nimport pandas as pd\nfrom numba import njit\nfrom pycocotools import _mask as coco_mask\nfrom sklearn.metrics import average_precision_score\nfrom sklearn.model_selection import KFold\nfrom tqdm import tqdm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load and align the training data\ninput_df = pd.read_parquet(\"train_features.parquet\").set_index(\"key\")\ntarget_df = pd.read_csv(\"../input/hpa-single-cell-image-classification/train.csv\")\nexpanded_labels = target_df[\"Label\"].str.split(\"|\", expand=True).astype(\"float\")\ntarget_classes = pd.Series(reduce(np.union1d, expanded_labels.apply(pd.unique))).dropna().astype(\"int\")\nmultihot_labels = pd.DataFrame(\n    {\n        tc: expanded_labels.eq(tc).any(axis=1)\n        for tc in target_classes\n    }\n)\ntarget_df = pd.concat([target_df[\"ID\"].rename(\"key\"), multihot_labels], axis=1).set_index(\"key\")\ntarget_df = target_df.loc[input_df.index]\nprint(\"Input data excerpt\")\ndisplay(input_df.iloc[:3, :20].T)\nprint(\"Target data excerpt\")\ndisplay(target_df.head(3).T)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load test data\ntest_input_df = pd.read_parquet(\"test_features.parquet\").set_index(\"key\")\nprint(\"Test input excerpt\")\ndisplay(test_input_df.iloc[:3, :20].T)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load and re-encode masks\n\ndef rle_to_mask(mask_nrow, mask_ncol, run_values, run_lengths):\n    \"\"\"NOTE: even with itertools, this is super slow and takes ~1sec per mask.\"\"\"\n    return np.array(tuple(itertools.chain.from_iterable(\n        itertools.repeat(v, r)\n        for v, r in zip(run_values, run_lengths)\n    ))).reshape((mask_nrow, mask_ncol), order=\"F\")\n\n\n@njit\ndef _rle_to_mask_1d(length, run_values, run_lengths):\n    \"\"\"Python loops are slow, but LLVM-compiled with numba is fast.\n    (~100x faster than `rle_to_mask`)\n    \"\"\"\n    res = np.empty(length, dtype=np.uint8)\n    i = 0\n    for v, r in zip(run_values, run_lengths):\n        for _ in range(r):\n            res[i] = v\n            i += 1\n    return res\n\n\ndef fast_rle_to_mask(mask_nrow, mask_ncol, run_values, run_lengths):\n    # note: it's faster to convert the lists to numpy arrays\n    #     rather than having the numba-jitted code take in Python lists\n    length_1d = mask_nrow * mask_ncol\n    res = _rle_to_mask_1d(length_1d, np.array(run_values), np.array(run_lengths))\n    return res.reshape((mask_nrow, mask_ncol), order=\"F\")\n\n\ndef binary_mask_to_ascii(mask):\n    # convert input mask to expected COCO API input --\n    mask_to_encode = mask.reshape(mask.shape[0], mask.shape[1], 1)\n    mask_to_encode = mask_to_encode.astype(np.uint8)\n    mask_to_encode = np.asfortranarray(mask_to_encode)\n\n    # RLE encode mask --\n    encoded_mask = coco_mask.encode(mask_to_encode)[0][\"counts\"]\n\n    # compress and base64 encoding --\n    binary_str = zlib.compress(encoded_mask, zlib.Z_BEST_COMPRESSION)\n    base64_str = base64.b64encode(binary_str)\n    return base64_str.decode()\n\n\ndef reencode_mask(mask_dict):\n    return binary_mask_to_ascii(fast_rle_to_mask(**mask_dict))\n\n\ncompetition_mask_strings = []\nimage_width_heights = []\nwith multiprocessing.Pool(multiprocessing.cpu_count()) as pool:\n    for mask_file_path in sorted(Path(\"rle_masks/\").iterdir()):\n        print(f\"Parsing mask file {mask_file_path}\")\n        with mask_file_path.open() as f:\n            rle_mask_json = json.load(f)\n        reencoded_masks = pool.imap(reencode_mask, tqdm(rle_mask_json))\n        competition_mask_strings.extend(reencoded_masks)\n        image_width_heights.extend([(mask[\"mask_ncol\"], mask[\"mask_nrow\"]) for mask in rle_mask_json])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Train one model per class\n\n# # big model for accuracy\n# model_params = dict(\n#     n_estimators=1_000,\n#     learning_rate=0.02,\n#     num_leaves=127,\n#     colsample_bytree=0.1,\n#     subsample=0.9,\n#     subsample_freq=1,\n# )\n\n# small model for testing\nmodel_params = dict(\n    n_estimators=30,\n    learning_rate=0.3,\n    num_leaves=7,\n    colsample_bytree=0.1,\n    subsample=0.9,\n    subsample_freq=1,\n)\n\n# train models on each class\nX = input_df.to_numpy()\nmodels = []\nfor i, target_class in enumerate(tqdm(target_df.columns)):\n    y = target_df[target_class].to_numpy()\n    model = lgb.LGBMClassifier(**model_params)\n    model.fit(X, y)\n    models.append(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# run predictions\nX = test_input_df.to_numpy()\npred_df = pd.DataFrame(\n    {\n        target_class: model.predict_proba(X)[:, 1]\n        for target_class, model in zip(target_df.columns, models)\n    },\n    index=test_input_df.index\n)\nprint(\"Prediction excerpt\")\npred_df.head(3).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create submission\n\n# let's choose thresholds so that the rate we guess a particular label is\n# a fixed multiple of the class frequency\nthresholds = {}\nrate_multiple = 2\nbase_rates = target_df.mean()\nthreshold_rates = rate_multiple * base_rates\nfor target_class, predictions in pred_df.iteritems():\n    quantile = min(1, threshold_rates[target_class])\n    thresholds[target_class] = predictions.quantile(quantile)\ndisplay(pd.DataFrame(dict(\n    base_rates=base_rates, threshold_rates=threshold_rates, threshold_values=thresholds\n)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_rows = []\nfor mask_string, (width, height), pred_tuple in zip(\n    tqdm(competition_mask_strings),\n    image_width_heights,\n    pred_df.itertuples(index=True, name=None)\n):\n    key, *pred = pred_tuple\n    for target_class, pred in zip(target_df.columns, pred):\n        if pred > thresholds[target_class]:\n            prediction_string = f\"{target_class} {pred:.20f} {mask_string}\"\n            row = {\n                \"ID\": key,\n                \"ImageWidth\": width,\n                \"ImageHeight\": height,\n                \"PredictionString\": prediction_string\n                \n            }\n            submission_rows.append(row)\nsubmission_df = pd.DataFrame.from_records(submission_rows)\nsubmission_df.to_csv(\"submission.csv\", index=False)\nsubmission_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Conclusion\n\nSo here we have it, no pre-training, no GPU, all thanks to Julia.\n\nI hope this work is inspiring to the Julia community. I've been a career-long Python programmer, and this is only my second project with the language (my first being a [library to download and parse stock data from IEX](https://github.com/lukemerrick/InvestorsExchange.jl)). Already I feel that I can try out so many things that just aren't feasible in Python/numpy, like non-ML segmentation algorithms, hand-implemented logistic regression code, and building custom visualizations of images. To any more veteran members of the Julia community who are reading through, I would deeply appreciate any feedback and constructive criticism on my code. I'm still a beginner trying to learn, and every tip helps!\n\nI also hope this work is inspiring to folks who aren't afraid to think creatively in the face of cookie-cutter deep learning model frameworks. Certainly I doubt my segmentation code gives results as clean as those from the offical [HPA Cell Segmentation](https://github.com/CellProfiling/HPA-Cell-Segmentation) project, but it gives solid results and runs quite efficiently (<1sec to segment a 400x400-resized image on a single core of a laptop CPU, most of the runtime is actually from doing higher-resolution feature computation). Not to devalue all the creativity and endless tweaking that goes into making a competitive deep learning submission, but I hope more folks out there will try more \"way out there\" approaches like this.\n\n\n### Opportunities for future work\n\nCollaboration welcome!\n\n- The first thing I want to do is go back and explain what all the Julia code does, with lots of images showing the intermediate steps of segmentation, spatial feature engineering, etc.\n\n- When it comes to improving my submission, I think there are some fairly simple things to try adding to the GBM modeling. \n    - It would be good to take into account the correlation between classes and the fact taht when multiple classes are present together, the features look different than when the classes are present individually. One thing to try would be finding the top target interactions (i.e. pairs of classes that show up together) and training a model that predicts the presence/absence of the interaction.\n    - Filtering the training data would be good, too, since a lot of cells in images may not actually match the label for the whole image. One simple approach would be to get a set of cross-validated predictions using the full data, and then drop the labels from cells whcih are predicted with significantly lower confidence than the other cells in the same image.\n\n- A less unconventional next step would be to get some autograd and convolutions going, but in Julia rather than Python. [Flux.jl](https://github.com/FluxML/Flux.jl) seems like a fun time, but also like another big chunk of work to dive into, so I'll have to see.\n\n- There's also the interesting technical challenge of seeing how much of my Julia code I can port over to run on GPU. It would probably require re-implementing the underlying Watershed algorithm from [ImageSegmentation.jl](https://juliaimages.org/v0.20/imagesegmentation/), which seems both daunting and interesting. It seems [possible](https://github.com/louismullie/watershed-cuda)."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}