*! version 1.2.3 20261030 github.com/mdroste/stata-ctools

program define cbinscatter, eclass
    version 14.1

    * Check observation limit (Stata plugin API limitation)
    if _N > 2147483647 {
        di as error "ctools does not support datasets exceeding 2^31 (2.147 billion) observations"
        di as error "This is a limitation of Stata's plugin API"
        exit 920
    }

    syntax varlist(min=2 max=2 fv) [aw fw pw iw] [if] [in], ///
        [NQuantiles(integer 20)] ///
        [Controls(varlist fv)] ///
        [Absorb(varlist)] ///
        [BY(varname)] ///
        [LINEtype(string)] ///
        [DISCrete] ///
        [GENxq(name)] ///
        [SAVEdata(string)] ///
        [REPORTreg] ///
        [NOgraph] ///
        [Verbose] ///
        [METHod(string)] ///
        [THReads(integer 0)] ///
        /* Graph options */ ///
        [TITLE(string)] ///
        [YTItle(string)] ///
        [XTItle(string)] ///
        [LEGEND(string asis)] ///
        [COLors(string)] ///
        [MSYMbols(string)] ///
        [MLAbels(string asis)] ///
        [*]

    * =========================================================================
    * UPFRONT VALIDATION - check all options before any data manipulation
    * This ensures errors don't leave the dataset in a bad state
    * =========================================================================

    * Own only the frame created by this invocation
    tempname graph_frame

    * Validate nquantiles early
    if `nquantiles' < 0 {
        di as error "nquantiles() invalid"
        exit 3300
    }
    if `nquantiles' < 2 | `nquantiles' > 1000 {
        di as error "cbinscatter: nquantiles must be between 2 and 1000"
        exit 198
    }

    * Validate linetype option early
    local linetype_num = 1   // Default to linear
    if "`linetype'" == "" {
        local linetype = "linear"
        local linetype_lower = "linear"
        local linetype_num = 1
    }
    else {
        local linetype_lower = lower("`linetype'")
        if "`linetype_lower'" == "none" | "`linetype_lower'" == "n" {
            local linetype_num = 0
        }
        else if "`linetype_lower'" == "linear" | "`linetype_lower'" == "line" | "`linetype_lower'" == "lfit" {
            local linetype_num = 1
        }
        else if "`linetype_lower'" == "quadratic" | "`linetype_lower'" == "qfit" {
            local linetype_num = 2
        }
        else if "`linetype_lower'" == "cubic" {
            local linetype_num = 3
        }
        else if "`linetype_lower'" == "connect" {
            local linetype_num = 0  // No fit, just connected points
        }
        else {
            di as error "cbinscatter: invalid linetype() option: `linetype'"
            di as error "  valid options: none, linear, qfit, cubic, connect"
            exit 198
        }
    }

    if "`reportreg'" != "" & `linetype_num' == 0 {
        di as error "Cannot specify reportreg when no fit line is being created."
        exit 198
    }
    local connect_opt
    if "`linetype_lower'" == "connect" local connect_opt "connect(l)"

    * Validate method option early (Cattaneo et al. "On Binscatter")
    * 0 = classic (default): residualize both Y and X, bin on residualized X
    * 1 = binsreg: bin on raw X, residualize Y only (conditional means)
    local method_num = 0   // Default to classic
    if "`method'" != "" {
        local method_lower = lower("`method'")
        if "`method_lower'" == "classic" | "`method_lower'" == "c" | "`method_lower'" == "residualize" {
            local method_num = 0
        }
        else if "`method_lower'" == "binsreg" | "`method_lower'" == "b" | "`method_lower'" == "conditional" {
            local method_num = 1
        }
        else {
            di as error "cbinscatter: invalid method() option: `method'"
            di as error "  valid options: classic, binsreg"
            exit 198
        }
    }

    * Validate savedata filename if specified (check directory exists)
    if "`savedata'" != "" {
        * Extract directory path
        local savedir = ""
        local lastslash = max(strpos(reverse("`savedata'"), "/"), strpos(reverse("`savedata'"), "\"))
        if `lastslash' > 0 {
            local savedir = substr("`savedata'", 1, strlen("`savedata'") - `lastslash')
        }
        if "`savedir'" != "" {
            capture confirm file "`savedir'/."
            if _rc {
                di as error "cbinscatter: directory does not exist for savedata(): `savedir'"
                exit 601
            }
        }
    }

    if "`genxq'" != "" {
        di as error "cbinscatter: genxq() is not implemented"
        exit 198
    }

    * =========================================================================
    * END UPFRONT VALIDATION - now proceed with data processing
    * =========================================================================

    local weight_var ""
    local weight_type = 0

    * Parse variable list
    gettoken depvar xvar : varlist
    local depvar = trim("`depvar'")
    local xvar = trim("`xvar'")

    * Validate that y and x are numeric (not string) - match binscatter rc=109
    capture confirm numeric variable `depvar'
    if _rc {
        di as error "type mismatch"
        exit 109
    }
    capture confirm numeric variable `xvar'
    if _rc {
        di as error "type mismatch"
        exit 109
    }

    local num_absorb : word count `absorb'

    * Count the requested sample before excluding incomplete observations
    marksample touse, novarlist

    * Stata counts the sample and drops incomplete observations only where
    * the ado needs them before the plugin runs (by-groups, discrete bins,
    * the verbose header); otherwise the plugin does both while it loads
    local stata_sample = ("`by'" != "" | "`discrete'" != "" | "`verbose'" != "")
    if `stata_sample' {
        qui count if `touse'
        local nobs = r(N)
        if `nobs' == 0 {
            di as error "no observations"
            exit 3301
        }
    }

    local check_weights = 0
    if "`weight'" != "" {
        * A weight that is a numeric variable is passed as it is, and the
        * plugin checks it; other expressions are evaluated and checked here
        local wexp = strtrim(substr(`"`exp'"', 2, .))
        capture unab wvar : `wexp', max(1)
        if !_rc capture confirm numeric variable `wvar'
        if !_rc & "`wvar'" == "`wexp'" {
            local weight_var "`wvar'"
            local check_weights = 1
        }
        else {
            tempvar evaluated_weight
            _ctools_weight, generate(`evaluated_weight') touse(`touse') type(`weight') expression(`"`exp'"')
            local weight_var "`evaluated_weight'"
        }
        if "`weight'" == "aweight" local weight_type = 1
        else if "`weight'" == "fweight" local weight_type = 2
        else if "`weight'" == "pweight" local weight_type = 3
        else if "`weight'" == "iweight" local weight_type = 4
    }
    local controls_original `controls'
    if "`controls'" != "" {
        fvrevar `controls' if `touse'
        local controls `r(varlist)'
    }
    local num_controls : word count `controls'

    * Group IDs must be dense after excluding incomplete observations.
    * absorb() values are labels (numeric or string); the plugin maps each
    * variable to dense levels over the final sample.
    if `stata_sample' {
        markout `touse' `depvar' `xvar' `controls' `weight_var'
        if "`absorb'" != "" markout `touse' `absorb', strok
        if "`by'" != "" markout `touse' `by', strok
    }
    local num_by_groups = 1
    if "`by'" != "" {
        local by_orig "`by'"
        tempvar by_numeric
        quietly _ctools_group `by' if `touse', generate(`by_numeric') label
        local by "`by_numeric'"
        local num_by_groups = r(N_groups)
        if `num_by_groups' == 0 error 2000
    }

    * Load the platform-appropriate ctools plugin if not already loaded
    _ctools_load
    * Stata scopes plugin registrations to the calling ado program.
    capture program ctools_plugin, plugin using("`__ctools_plugin'")
    if _rc != 0 & _rc != 110 exit 601
    capture confirm number 0

    * Set up parameters via Stata scalars
    scalar __cbinscatter_nquantiles = `nquantiles'
    scalar __cbinscatter_linetype = `linetype_num'
    scalar __cbinscatter_compute_se = 0
    scalar __cbinscatter_discrete = ("`discrete'" != "")
    scalar __cbinscatter_has_controls = (`num_controls' > 0)
    scalar __cbinscatter_num_controls = `num_controls'
    scalar __cbinscatter_has_absorb = (`num_absorb' > 0)
    scalar __cbinscatter_num_absorb = `num_absorb'
    scalar __cbinscatter_has_by = ("`by'" != "")
    scalar __cbinscatter_has_weights = (`weight_type' > 0)
    scalar __cbinscatter_weight_type = `weight_type'
    scalar __cbinscatter_verbose = ("`verbose'" != "")
    scalar __cbinscatter_reportreg = ("`reportreg'" != "")
    scalar __cbinscatter_method = `method_num'
    scalar __cbinscatter_maxiter = 10000
    scalar __cbinscatter_tolerance = 1e-10
    * The plugin selects the sample from a marker variable (or, without
    * if/in, from the data alone) rather than from an if condition, which
    * Stata answers one observation at a time
    scalar __cbinscatter_touse_mode = cond(`"`if'`in'"' != "", 1, 2)
    scalar __cbinscatter_check_weights = `check_weights'

    * Display info if verbose
    if "`verbose'" != "" {
        di as text ""
        di as text "{hline 60}"
        di as text "cbinscatter: C-Accelerated Binned Scatter"
        di as text "{hline 60}"
        di as text "Y variable: " as result "`depvar'"
        di as text "X variable: " as result "`xvar'"
        if `num_controls' > 0 {
            di as text "Controls:   " as result "`controls'"
        }
        if `num_absorb' > 0 {
            di as text "Absorb:     " as result "`absorb'"
        }
        if "`by'" != "" {
            di as text "By:         " as result "`by'"
        }
        di as text "Bins:       " as result `nquantiles'
        di as text "Obs:        " as result `nobs'
        di as text "{hline 60}"
        di ""
    }

    * Create matrices for C plugin to fill
    local max_bins = `nquantiles'
    if "`discrete'" != "" {
        * For discrete, could have many more bins - use a reasonable max
        local max_bins = min(`nobs', 500)
    }
    matrix __cbinscatter_bins = J(`max_bins' * `num_by_groups', 5, .)
    if `linetype_num' > 0 {
        matrix __cbinscatter_coefs = J(`num_by_groups', 4, .)
        matrix __cbinscatter_fit_stats = J(`num_by_groups', 2, .)
    }

    * Build varlist for plugin: y x [controls] [absorb] [by] [weight]
    local plugin_varlist `depvar' `xvar'
    if "`controls'" != "" {
        local plugin_varlist `plugin_varlist' `controls'
    }
    if "`absorb'" != "" {
        local plugin_varlist `plugin_varlist' `absorb'
    }
    if "`by'" != "" {
        local plugin_varlist `plugin_varlist' `by'
    }
    if `weight_type' > 0 {
        local plugin_varlist `plugin_varlist' `weight_var'
    }
    if `"`if'`in'"' != "" {
        local plugin_varlist `plugin_varlist' `touse'
    }

    * Build threads option string
    local threads_code ""
    if `threads' > 0 {
        local threads_code "threads(`threads')"
    }

    * Record start time
    quietly _ctools_clock __ct_lap99

    * Call the C plugin
    capture noisily plugin call ctools_plugin `plugin_varlist', ///
        "cbinscatter `threads_code' compute_bins"

    local plugin_rc = _rc
    if `plugin_rc' {
        * The plugin itself reports an empty sample and invalid weights
        capture local reported = __cbinscatter_reported
        if _rc local reported = 0
        _ctools_cleanup cbinscatter
        if !`reported' di as error "cbinscatter: plugin error (rc=`plugin_rc')"
        exit `plugin_rc'
    }
    if !`stata_sample' local nobs = __cbinscatter_N_touse

    _ctools_lap elapsed __ct_lap99 `__ct_lap99'

    * Retrieve results from scalars
    local N_used = __cbinscatter_N
    local N_dropped = `nobs' - `N_used'
    local actual_num_groups = __cbinscatter_num_groups

    * reportreg prints the underlying fit regressions, as binscatter does.
    * Compute these only on request; plotted estimates remain the plugin's.
    if "`reportreg'" != "" {
        local ry `depvar'
        local rx `xvar'
        local wt
        if "`weight'" != "" local wt [`weight'=`weight_var']
        * areg retains singleton FE groups, matching binscatter's residualization.
        * Encode labels for string/negative FE values; additional FEs are dummies.
        local report_absorb
        local report_fe
        foreach fe of local absorb {
            tempvar report_group
            quietly _ctools_group `fe' if `touse', generate(`report_group')
            if "`report_absorb'" == "" local report_absorb `report_group'
            else local report_fe `report_fe' i.`report_group'
        }
        if `method_num' == 0 & "`controls'`absorb'" != "" {
            foreach role in ry rx {
                tempvar resid
                if "`absorb'" == "" {
                    quietly regress ``role'' `controls' `wt' if `touse'
                    quietly predict double `resid' if `touse', residuals
                }
                else {
                    quietly areg ``role'' `controls' `report_fe' `wt' if `touse', absorb(`report_absorb')
                    quietly predict double `resid' if `touse', residuals
                }
                quietly summarize ``role'' `wt' if `touse', meanonly
                quietly replace `resid' = `resid' + r(mean) if `touse'
                local `role' `resid'
            }
        }
        local fitvars `rx'
        if `linetype_num' >= 2 {
            tempvar quadratic
            quietly gen double `quadratic' = `rx'^2 if `touse'
            local fitvars `quadratic' `fitvars'
        }
        if `linetype_num' == 3 {
            tempvar cubic
            quietly gen double `cubic' = `rx'^3 if `touse'
            local fitvars `cubic' `fitvars'
        }
        forvalues g=1/`actual_num_groups' {
            local group_if
            if "`by'" != "" local group_if & `by' == `g'
            di as text "Fit regression: `depvar' on `xvar' (group `g')"
            if `method_num' == 1 & "`absorb'" != "" {
                areg `ry' `fitvars' `controls' `report_fe' `wt' if `touse' `group_if', absorb(`report_absorb')
            }
            else {
                local fitcontrols
                if `method_num' == 1 local fitcontrols `controls'
                regress `ry' `fitvars' `fitcontrols' `wt' if `touse' `group_if'
            }
        }
    }

    * Display timing if requested
    if "`verbose'" != "" {
        local plugin_overhead = `elapsed' - __cbinscatter_time_total
        if `plugin_overhead' < 0 local plugin_overhead = 0
        di as text ""
        di as text "{hline 55}"
        di as text "cbinscatter timing breakdown:"
        di as text "{hline 55}"
        di as text "  C plugin internals:"
        di as text "    Data load:              " as result %8.4f __cbinscatter_time_load " sec"
        di as text "    Residualize:            " as result %8.4f __cbinscatter_time_resid " sec"
        di as text "    Bin computation:        " as result %8.4f __cbinscatter_time_bins " sec"
        di as text "    Line fitting:           " as result %8.4f __cbinscatter_time_fit " sec"
        di as text "  {hline 53}"
        di as text "    C plugin total:         " as result %8.4f __cbinscatter_time_total " sec"
        di as text "  {hline 53}"
        di as text "  Plugin call overhead:     " as result %8.4f `plugin_overhead' " sec"
        di as text "{hline 55}"
        di as text "    Wall clock total:       " as result %8.4f `elapsed' " sec"
        di as text "{hline 55}"

        * Display thread diagnostics
        capture local __threads_max = __cbinscatter_threads_max
        if _rc == 0 {
            capture local __openmp_enabled = __cbinscatter_openmp_enabled
            if _rc != 0 local __openmp_enabled = 0
            di as text ""
            di as text "  Thread diagnostics:"
            di as text "    OpenMP enabled:         " as result %8.0f `__openmp_enabled'
            di as text "    Max threads available:  " as result %8.0f `__threads_max'
            di as text "{hline 55}"
        }
    }

    * Display summary
    if `N_dropped' > 0 {
        di as text "(dropped " as result `N_dropped' as text " observations due to missing values)"
    }

    * Save bin data if requested. The bins are built in a frame of their own
    * (Stata 16 or later), so the data in memory are neither copied nor
    * cleared; before Stata 16 they are preserved (on disk above
    * c(max_preservemem)) and restored.
    if "`savedata'" != "" {
        if c(stata_version) >= 16 {
            tempname binframe
            frame create `binframe'
            frame `binframe': _cbinscatter_savebins `"`savedata'"'
            frame drop `binframe'
        }
        else {
            preserve
            _cbinscatter_savebins `"`savedata'"'
            restore
        }
        di as text "(bin data saved to " as result "`savedata'" as text ")"
    }

    * Generate bin assignment variable if requested


    * Create graph if not suppressed
    if "`nograph'" == "" {
        * Default colors (same as binscatter)
        local colorlist "navy maroon forest_green dkorange teal cranberry lavender khaki sienna emidblue emerald brown erose gold bluishgray"

        * Parse user color options or use defaults
        if "`colors'" == "" {
            local colors "`colorlist'"
        }

        * Set marker colors and line colors (binscatter uses separate defaults)
        * For single y-var with single by-group (no connect): mcolors=color1, lcolors=color2
        * Otherwise: both use the same color sequence
        local mcolors "`colors'"
        local lcolors "`colors'"
        if `actual_num_groups' == 1 & "`linetype_lower'" != "connect" {
            local mcolors : word 1 of `colors'
            local lcolors : word 2 of `colors'
        }

        * Get bin data from matrix
        tempname binmat
        matrix `binmat' = __cbinscatter_bins

        * Find x range for fit lines
        local x_min = .
        local x_max = .
        local nrows = rowsof(`binmat')
        forval r = 1/`nrows' {
            local xval = `binmat'[`r', 3]
            if `xval' != . {
                if `x_min' == . | `xval' < `x_min' {
                    local x_min = `xval'
                }
                if `x_max' == . | `xval' > `x_max' {
                    local x_max = `xval'
                }
            }
        }

        * Build scatter commands using scatteri (like binscatter)
        local scatters ""
        local fits ""
        local legend_labels ""
        local legend_order ""

        if `actual_num_groups' == 1 {
            * Single group - build scatteri coordinate list
            local coords ""
            * scatteri labels each point with text that follows its coordinates
            local marker_label ""
            if `"`mlabels'"' != "" {
                gettoken marker_label : mlabels
                local marker_label `""`marker_label'""'
            }
            forval r = 1/`nrows' {
                local yval = `binmat'[`r', 4]
                local xval = `binmat'[`r', 3]
                if `yval' != . & `xval' != . {
                    local coords `"`coords' `yval' `xval' `marker_label'"'
                }
            }

            * Get marker and line colors (binscatter uses different colors for single series)
            local thismcolor : word 1 of `mcolors'
            local thislcolor : word 1 of `lcolors'
            local thissymbol ""
            local symbol : word 1 of `msymbols'
            if "`symbol'" != "" local thissymbol "msymbol(`symbol')"

            * Build scatter with mcolor and lcolor (lcolor for connect mode)
            local scatters `"(scatteri `coords', mcolor(`thismcolor') lcolor(`thislcolor') `thissymbol' `connect_opt')"'

            * Add fit line if requested
            if `linetype_num' > 0 {
                local coef_const = __cbinscatter_coefs[1, 1]
                local coef_linear = __cbinscatter_coefs[1, 2]

                if `linetype_num' == 1 {
                    * Linear fit using function
                    local fits "(function `coef_linear'*x+`coef_const', range(`x_min' `x_max') lcolor(`thislcolor'))"
                }
                else if `linetype_num' >= 2 {
                    local coef_quad = __cbinscatter_coefs[1, 3]
                    if `linetype_num' == 2 {
                        local fits "(function `coef_quad'*x^2+`coef_linear'*x+`coef_const', range(`x_min' `x_max') lcolor(`thislcolor'))"
                    }
                    else {
                        local coef_cubic = __cbinscatter_coefs[1, 4]
                        local fits "(function `coef_cubic'*x^3+`coef_quad'*x^2+`coef_linear'*x+`coef_const', range(`x_min' `x_max') lcolor(`thislcolor'))"
                    }
                }
            }

            * Legend off for single group
            local legend_opt "legend(off)"
        }
        else {
            * Multiple by-groups
            local legend_idx = 1

            forval g = 1/`actual_num_groups' {
                * Build coordinate list for this group
                local coords ""
                local grp_x_min = .
                local grp_x_max = .
                local marker_label ""
                if `"`mlabels'"' != "" {
                    local marker_label : word `g' of `mlabels'
                    if `"`marker_label'"' != "" local marker_label `""`marker_label'""'
                }

                forval r = 1/`nrows' {
                    local grpval = `binmat'[`r', 1]
                    if `grpval' == `g' {
                        local yval = `binmat'[`r', 4]
                        local xval = `binmat'[`r', 3]
                        if `yval' != . & `xval' != . {
                            local coords `"`coords' `yval' `xval' `marker_label'"'
                            if `grp_x_min' == . | `xval' < `grp_x_min' {
                                local grp_x_min = `xval'
                            }
                            if `grp_x_max' == . | `xval' > `grp_x_max' {
                                local grp_x_max = `xval'
                            }
                        }
                    }
                }

                * Get marker and line colors for this series
                local thismcolor : word `g' of `mcolors'
                local thislcolor : word `g' of `lcolors'
                if "`thismcolor'" == "" {
                    local thismcolor "navy"
                }
                if "`thislcolor'" == "" {
                    local thislcolor "navy"
                }

                local thissymbol ""
                local symbol : word `g' of `msymbols'
                if "`symbol'" != "" local thissymbol "msymbol(`symbol')"

                local scatters `"`scatters' (scatteri `coords', mcolor(`thismcolor') lcolor(`thislcolor') `thissymbol' `connect_opt')"'
                local group_label : label (`by') `g'
                local legend_labels `"`legend_labels' label(`legend_idx' "`by_orig'==`group_label'")"'
                local legend_order "`legend_order' `legend_idx'"
                local legend_idx = `legend_idx' + 1

                * Add fit line for this group if requested
                if `linetype_num' > 0 & `"`coords'"' != "" {
                    local coef_const = __cbinscatter_coefs[`g', 1]
                    local coef_linear = __cbinscatter_coefs[`g', 2]

                    if `linetype_num' == 1 {
                        local fits "`fits' (function `coef_linear'*x+`coef_const', range(`grp_x_min' `grp_x_max') lcolor(`thislcolor'))"
                    }
                    else if `linetype_num' >= 2 {
                        local coef_quad = __cbinscatter_coefs[`g', 3]
                        if `linetype_num' == 2 {
                            local fits "`fits' (function `coef_quad'*x^2+`coef_linear'*x+`coef_const', range(`grp_x_min' `grp_x_max') lcolor(`thislcolor'))"
                        }
                        else {
                            local coef_cubic = __cbinscatter_coefs[`g', 4]
                            local fits "`fits' (function `coef_cubic'*x^3+`coef_quad'*x^2+`coef_linear'*x+`coef_const', range(`grp_x_min' `grp_x_max') lcolor(`thislcolor'))"
                        }
                    }
                }
            }

            * Build legend option
            if `"`legend'"' == "" {
                local legend_opt `"legend(`legend_labels' order(`legend_order'))"'
            }
            else {
                local legend_opt `"legend(`legend')"'
            }
        }

        * Set default titles if not specified (like binscatter)
        if `"`ytitle'"' == "" {
            local ytitle "`depvar'"
        }
        if `"`xtitle'"' == "" {
            local xtitle "`xvar'"
        }

        * Build title option only if specified
        local title_opt ""
        if `"`title'"' != "" {
            local title_opt `"title(`"`title'"')"'
        }

        * Execute graph (matching binscatter format)
        * Use temporary frame to avoid twoway overhead on large datasets
        * (twoway speed depends on dataset size in memory, even for scatteri)
        local _use_frame = 0
        if c(stata_version) >= 16 {
            local _use_frame = 1
        }

        if `_use_frame' {
            * Stata 16+: use frames for speed (no preserve/restore overhead)
            * Use capture to ensure proper cleanup if twoway fails
            local _orig_frame = c(frame)
            capture frame drop `graph_frame'
            frame create `graph_frame'
            frame change `graph_frame'

            capture noisily twoway `scatters' `fits', ///
                graphregion(fcolor(white)) ///
                ytitle(`"`ytitle'"') ///
                xtitle(`"`xtitle'"') ///
                `legend_opt' ///
                `title_opt' ///
                `options'

            local _graph_rc = _rc

            * Always clean up frame, even on error
            frame change `_orig_frame'
            capture frame drop `graph_frame'

            if `_graph_rc' {
                di as error "cbinscatter: graph generation failed (rc=`_graph_rc')"
                exit `_graph_rc'
            }
        }
        else {
            * Stata 14-15: fall back to preserve/keep/restore
            * Use capture to ensure restore happens even on error
            preserve
            qui keep in 1/1

            capture noisily twoway `scatters' `fits', ///
                graphregion(fcolor(white)) ///
                ytitle(`"`ytitle'"') ///
                xtitle(`"`xtitle'"') ///
                `legend_opt' ///
                `title_opt' ///
                `options'

            local _graph_rc = _rc
            restore

            if `_graph_rc' {
                di as error "cbinscatter: graph generation failed (rc=`_graph_rc')"
                exit `_graph_rc'
            }
        }
    }

    * Post e() results
    ereturn clear
    ereturn scalar N = `N_used'
    ereturn scalar N_dropped = `N_dropped'
    ereturn scalar nquantiles = `nquantiles'
    ereturn scalar num_groups = `actual_num_groups'
    ereturn matrix bindata = __cbinscatter_bins
    if `linetype_num' > 0 {
        ereturn matrix coefs = __cbinscatter_coefs
        ereturn matrix fit_stats = __cbinscatter_fit_stats
    }
    ereturn local depvar "`depvar'"
    ereturn local xvar "`xvar'"
    if "`controls'" != "" {
        ereturn local controls "`controls_original'"
    }
    if "`absorb'" != "" {
        ereturn local absorb "`absorb'"
    }
    if "`by'" != "" {
        ereturn local by "`by_orig'"
    }
    ereturn local linetype "`linetype'"
    ereturn local cmd "cbinscatter"
    ereturn local cmdline "cbinscatter `0'"

    _ctools_cleanup cbinscatter
end

* Replace the current dataset with the bins in matrix __cbinscatter_bins and
* save them.
program define _cbinscatter_savebins
    version 14.1
    args savedata
    quietly {
        clear
        svmat __cbinscatter_bins, names(col)
        rename c1 by_group
        rename c2 bin_id
        rename c3 x_mean
        rename c4 y_mean
        rename c5 n_obs
        * Drop empty rows
        drop if missing(bin_id)
        save `"`savedata'"', replace
    }
end