#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <sstream>
#include <string>
#include <utility>
#include <vector>

struct FSegment
{
    int passage=0,index=0;
    std::array<double,3> a{},b{};
    double radiusA=0.0,radiusB=0.0;
};
struct FTile
{
    int strate=0,x=0,y=0,z=0,lod=0;
    std::uint64_t baselineFeatures=0,baselineCandidates=0,baselinePairs=0;
};
struct FTotals
{
    std::uint64_t tiles=0,featureTiles=0,proofTiles=0,groups=0,segments=0;
    std::uint64_t pairs=0,contributingPairs=0,touched=0,support=0;
    std::uint64_t overlapPairs=0,crossPassageOverlap=0;
    std::uint64_t baselineFeatureTiles=0,baselineFeatures=0,baselineCandidates=0;
    std::uint64_t baselinePairs=0;
};

static bool Intersects(const std::array<double,3>& amin,
    const std::array<double,3>& amax,const std::array<double,3>& bmin,
    const std::array<double,3>& bmax)
{
    for(int axis=0;axis<3;++axis)
        if(amax[axis]<bmin[axis]||amin[axis]>bmax[axis])return false;
    return true;
}

static int GridFirst(double bound,const double gridMin[3],double spacing,int axis)
{
    const auto boundSub=static_cast<std::int64_t>(std::floor(bound*1024.0));
    const auto gridMinSub=static_cast<std::int64_t>(std::llround(gridMin[axis]*1024.0));
    const auto stepSub=static_cast<std::int64_t>(std::llround(spacing*1024.0));
    const auto numerator=boundSub-gridMinSub;
    auto quotient=numerator/stepSub;
    if(numerator%stepSub>0)++quotient;
    return static_cast<int>(quotient);
}

static int GridLast(double bound,const double gridMin[3],double spacing,int axis)
{
    const auto boundSub=static_cast<std::int64_t>(std::ceil(bound*1024.0));
    const auto gridMinSub=static_cast<std::int64_t>(std::llround(gridMin[axis]*1024.0));
    const auto stepSub=static_cast<std::int64_t>(std::llround(spacing*1024.0));
    const auto numerator=boundSub-gridMinSub;
    auto quotient=numerator/stepSub;
    if(numerator%stepSub<0)--quotient;
    return static_cast<int>(quotient);
}

int main(int argc,char** argv)
{
    if(argc!=3)
    {
        std::cerr<<"usage: CountReachCandidates <input.tsv> <output.csv>\n";
        return 2;
    }
    std::ifstream input(argv[1]);
    if(!input){std::cerr<<"cannot open input\n";return 2;}
    std::vector<std::vector<FSegment>> passages(3);
    std::vector<FTile> tiles;
    std::string line;
    while(std::getline(input,line))
    {
        std::istringstream row(line);
        char kind=0;row>>kind;
        if(kind=='S')
        {
            FSegment s;row>>s.passage>>s.index
                >>s.a[0]>>s.a[1]>>s.a[2]>>s.b[0]>>s.b[1]>>s.b[2]
                >>s.radiusA>>s.radiusB;
            if(!row||s.passage<0||s.passage>=static_cast<int>(passages.size()))
            {std::cerr<<"bad segment row\n";return 2;}
            passages[static_cast<std::size_t>(s.passage)].push_back(s);
        }
        else if(kind=='T')
        {
            FTile t;row>>t.strate>>t.x>>t.y>>t.z>>t.lod
                >>t.baselineFeatures>>t.baselineCandidates>>t.baselinePairs;
            if(!row){std::cerr<<"bad tile row\n";return 2;}
            tiles.push_back(t);
        }
        else if(kind!='#'&&kind!='\0')
        {std::cerr<<"unknown row\n";return 2;}
    }
    for(auto& p:passages)
        std::sort(p.begin(),p.end(),[](const auto& a,const auto& b)
            {return a.index<b.index;});

    std::ofstream output(argv[2]);
    if(!output){std::cerr<<"cannot open output\n";return 2;}
    output<<"reach_m,bounds_mode,strate,lod,tiles,feature_tiles,proof_tiles,feature_groups,segment_pairs,contributing_pairs,touched_samples,round_support_samples,overlap_pairs,cross_passage_overlaps,baseline_feature_tiles,baseline_features,baseline_candidates,baseline_pairs\n";
    constexpr std::array<double,4> reaches={4.0,3.0,2.0,1.5};
    const std::array<std::pair<const char*,double>,2> boundsModes={{{"v225_core_voxel_units",0.25},
        {"physical_metres",1.0}}};
    constexpr double coreVoxelMetres=0.25;
    constexpr double marginMetres=0.01;
    for(const auto& [boundsMode,boundsScale]:boundsModes)
    for(const double reach:reaches)
    {
        const double requiredReach=std::max(reach+0.21+0.125,0.75+0.21);
        std::map<std::pair<int,int>,FTotals> totals;
        for(const auto& tile:tiles)
        {
            const int stepVoxels=1<<tile.lod;
            const double spacing=stepVoxels*coreVoxelMetres;
            const int dim=32+3;
            const double tileOrigin[3]={tile.x*32.0*spacing,tile.y*32.0*spacing,
                tile.z*32.0*spacing};
            const double gridMin[3]={tileOrigin[0]-spacing,tileOrigin[1]-spacing,
                tileOrigin[2]-spacing};
            const double gridMax[3]={gridMin[0]+(dim-1)*spacing,
                gridMin[1]+(dim-1)*spacing,gridMin[2]+(dim-1)*spacing};
            const std::array<double,3> queryMin={gridMin[0],gridMin[1],gridMin[2]};
            const std::array<double,3> queryMax={gridMax[0],gridMax[1],gridMax[2]};
            const std::array<double,3> proofMin={tileOrigin[0],tileOrigin[1],tileOrigin[2]};
            const std::array<double,3> proofMax={tileOrigin[0]+32.0*spacing,
                tileOrigin[1]+32.0*spacing,tileOrigin[2]+32.0*spacing};

            const std::size_t sampleCount=static_cast<std::size_t>(dim)*dim*dim;
            std::vector<std::vector<std::uint8_t>> touchedMasks(passages.size(),
                std::vector<std::uint8_t>(sampleCount,0));
            std::vector<std::vector<std::uint8_t>> supportMasks(passages.size(),
                std::vector<std::uint8_t>(sampleCount,0));
            std::vector<std::uint8_t> unionMask(sampleCount,0);
            FTotals local;
            bool anyFeature=false,anyProof=false;

            for(std::size_t p=0;p<passages.size();++p)
            {
                const auto& segments=passages[p];
                const std::size_t groupCount=(segments.size()+4)/5;
                for(std::size_t group=0;group<groupCount;++group)
                {
                    const std::size_t begin=group*5,end=std::min(begin+5,segments.size());
                    std::array<double,3> routeMin={
                        std::numeric_limits<double>::infinity(),
                        std::numeric_limits<double>::infinity(),
                        std::numeric_limits<double>::infinity()};
                    std::array<double,3> routeMax={
                        -std::numeric_limits<double>::infinity(),
                        -std::numeric_limits<double>::infinity(),
                        -std::numeric_limits<double>::infinity()};
                    double maxRadius=0.0;
                    for(std::size_t i=begin;i<end;++i)
                    {
                        const auto& s=segments[i];
                        maxRadius=std::max({maxRadius,s.radiusA,s.radiusB});
                        for(int axis=0;axis<3;++axis)
                        {
                            routeMin[axis]=std::min({routeMin[axis],s.a[axis],s.b[axis]});
                            routeMax[axis]=std::max({routeMax[axis],s.a[axis],s.b[axis]});
                        }
                    }
                    // The v225 implementation divides positions and radii by Q
                    // (Core voxels), then adds its *_Metres constants directly.
                    // Reproduce that code-unit behavior here for the baseline check.
                    const double featureExtent=maxRadius
                        +(requiredReach+marginMetres)*boundsScale;
                    for(int axis=0;axis<3;++axis)
                    {routeMin[axis]-=featureExtent;routeMax[axis]+=featureExtent;}
                    if(!Intersects(routeMin,routeMax,queryMin,queryMax))continue;
                    anyFeature=true;++local.groups;

                    for(std::size_t i=begin;i<end;++i)
                    {
                        const auto& s=segments[i];
                        const double extent=std::max(s.radiusA,s.radiusB)
                            +(requiredReach+marginMetres)*boundsScale;
                        std::array<double,3> low{},high{};
                        int first[3]{},last[3]{};
                        bool valid=true;
                        for(int axis=0;axis<3;++axis)
                        {
                            low[axis]=std::min(s.a[axis],s.b[axis])-extent;
                            high[axis]=std::max(s.a[axis],s.b[axis])+extent;
                            first[axis]=std::max(0,GridFirst(low[axis],gridMin,spacing,axis));
                            last[axis]=std::min(dim-1,GridLast(high[axis],gridMin,spacing,axis));
                            valid=valid&&first[axis]<=last[axis];
                        }
                        // Proof bounds are per segment and include the same conservative
                        // reach; they test the full tile box without the sample halo.
                        if(Intersects(low,high,proofMin,proofMax))anyProof=true;
                        if(!valid)continue;
                        ++local.segments;
                        for(int axis=0;axis<3;++axis)
                        {
                            const auto subLo=static_cast<std::int64_t>(std::floor(low[axis]*1024.0));
                            const auto subHi=static_cast<std::int64_t>(std::ceil(high[axis]*1024.0));
                            const auto gridLo=static_cast<std::int64_t>(std::llround(gridMin[axis]*1024.0));
                            const auto subStep=static_cast<std::int64_t>(std::llround(spacing*1024.0));
                            const auto expectedFirst=std::max<std::int64_t>(0,
                                (subLo-gridLo+subStep-1)/subStep);
                            const auto expectedLast=std::min<std::int64_t>(dim-1,
                                (subHi-gridLo)/subStep);
                            first[axis]=static_cast<int>(expectedFirst);
                            last[axis]=static_cast<int>(expectedLast);
                        }
                        if(first[0]>last[0]||first[1]>last[1]||first[2]>last[2])continue;
                        const std::uint64_t nx=static_cast<std::uint64_t>(last[0]-first[0]+1);
                        const std::uint64_t ny=static_cast<std::uint64_t>(last[1]-first[1]+1);
                        const std::uint64_t nz=static_cast<std::uint64_t>(last[2]-first[2]+1);
                        local.pairs+=nx*ny*nz;

                        for(int gz=first[2];gz<=last[2];++gz)
                        for(int gy=first[1];gy<=last[1];++gy)
                        for(int gx=first[0];gx<=last[0];++gx)
                        {
                            const double q[3]={gridMin[0]+gx*spacing,
                                gridMin[1]+gy*spacing,gridMin[2]+gz*spacing};
                            const double dx=s.b[0]-s.a[0],dy=s.b[1]-s.a[1],dz=s.b[2]-s.a[2];
                            const double lengthSq=dx*dx+dy*dy+dz*dz;
                            if(!(lengthSq>0.0))continue;
                            const double projection=(q[0]-s.a[0])*dx+(q[1]-s.a[1])*dy
                                +(q[2]-s.a[2])*dz;
                            const double t=std::clamp(projection/lengthSq,0.0,1.0);
                            const double ex=q[0]-(s.a[0]+t*dx);
                            const double ey=q[1]-(s.a[1]+t*dy);
                            const double ez=q[2]-(s.a[2]+t*dz);
                            const double radius=s.radiusA+t*(s.radiusB-s.radiusA);
                            const double distance=std::sqrt(ex*ex+ey*ey+ez*ez)-radius;
                            const double oneSubvoxel=1.0/1024.0;
                            if(distance<reach+oneSubvoxel)
                            {
                                ++local.contributingPairs;
                                const std::size_t index=(static_cast<std::size_t>(gz)*dim+gy)*dim+gx;
                                touchedMasks[p][index]=1;
                                unionMask[index]=1;
                                if(distance<=reach)supportMasks[p][index]=1;
                            }
                        }
                    }
                }
            }
            local.tiles=1;
            local.featureTiles=anyFeature?1:0;
            local.proofTiles=anyProof?1:0;
            local.baselineFeatureTiles=tile.baselineFeatures?1:0;
            local.baselineFeatures=tile.baselineFeatures;
            local.baselineCandidates=tile.baselineCandidates;
            local.baselinePairs=tile.baselinePairs;
            for(const auto& mask:touchedMasks)
                local.touched+=static_cast<std::uint64_t>(
                    std::count(mask.begin(),mask.end(),static_cast<std::uint8_t>(1)));
            for(std::size_t p=0;p<passages.size();++p)
            {
                const auto& mask=supportMasks[p];
                local.support+=static_cast<std::uint64_t>(
                    std::count(mask.begin(),mask.end(),static_cast<std::uint8_t>(1)));
            }
            const std::uint64_t unique=static_cast<std::uint64_t>(
                std::count(unionMask.begin(),unionMask.end(),static_cast<std::uint8_t>(1)));
            local.crossPassageOverlap=local.touched-unique;
            local.overlapPairs=local.contributingPairs-local.touched;
            auto& total=totals[{tile.strate,tile.lod}];
            total.tiles+=local.tiles;total.featureTiles+=local.featureTiles;
            total.proofTiles+=local.proofTiles;total.groups+=local.groups;
            total.segments+=local.segments;total.pairs+=local.pairs;
            total.contributingPairs+=local.contributingPairs;total.touched+=local.touched;
            total.support+=local.support;total.overlapPairs+=local.overlapPairs;
            total.crossPassageOverlap+=local.crossPassageOverlap;
            total.baselineFeatureTiles+=local.baselineFeatureTiles;
            total.baselineFeatures+=local.baselineFeatures;
            total.baselineCandidates+=local.baselineCandidates;
            total.baselinePairs+=local.baselinePairs;
        }
        for(const auto& [key,t]:totals)
            output<<std::fixed<<std::setprecision(2)<<reach<<','<<boundsMode<<','
                <<key.first<<','<<key.second
                <<','<<t.tiles<<','<<t.featureTiles<<','<<t.proofTiles<<','<<t.groups
                <<','<<t.pairs<<','<<t.contributingPairs<<','<<t.touched<<','<<t.support
                <<','<<t.overlapPairs<<','<<t.crossPassageOverlap<<','
                <<t.baselineFeatureTiles<<','<<t.baselineFeatures<<','
                <<t.baselineCandidates<<','<<t.baselinePairs<<'\n';
        std::cerr<<"[L129Precount] reach_m="<<reach<<" required_reach_m="<<requiredReach
            <<" grouped_rows="<<totals.size()<<"\n";
    }
    return 0;
}