# src/SciKit.jl
module SciKit

using LinearAlgebra
using Statistics

# 类型定义
include("types.jl")

# 子模块
include("utils.jl")
include("integration.jl")
include("differentiation.jl")
include("interpolation.jl")
include("optimization.jl")

# 导出公共 API
export
    # 积分
    quad, simpson, gauss_legendre,
    # 微分
    central_diff, gradient, jacobian,
    # 插值
    LinearInterp, CubicSpline,
    # 优化
    gradient_descent, newton_method,
    # 类型
    QuadratureResult

# 版本信息
const VERSION = v"0.1.0"

function __init__()
    @info "SciKit.jl v$VERSION loaded"
end

end # module

# src/types.jl

# 抽象类型层次
abstract type AbstractIntegrator end
abstract type AbstractDifferentiator end
abstract type AbstractInterpolator{T} end
abstract type AbstractOptimizer end

# 积分结果
"""
    QuadratureResult{T}

存储数值积分结果的结构体。

# 字段
- `value::T`: 积分值
- `error::T`: 估计误差
- `n_evals::Int`: 函数求值次数
- `converged::Bool`: 是否收敛
"""
struct QuadratureResult{T<:Real}
    value::T
    error::T
    n_evals::Int
    converged::Bool
end

# 参数化插值器
"""
    LinearInterp{T}

线性插值器，支持任意实数类型。
"""
struct LinearInterp{T<:Real} <: AbstractInterpolator{T}
    xs::Vector{T}
    ys::Vector{T}
    n::Int
    
    function LinearInterp(xs::Vector{T}, ys::Vector{T}) where T
        n = length(xs)
        length(ys) == n || throw(DimensionMismatch("xs 和 ys 长度不匹配"))
        issorted(xs) || throw(ArgumentError("xs 必须排序"))
        new{T}(xs, ys, n)
    end
end

"""
    CubicSpline{T}

三次样条插值器。
"""
struct CubicSpline{T<:Real} <: AbstractInterpolator{T}
    xs::Vector{T}
    ys::Vector{T}
    coeffs::Matrix{T}  # [a, b, c, d] 系数
    n::Int
end

# 优化结果
struct OptimizationResult{T}
    x_min::Vector{T}
    f_min::T
    n_iter::Int
    converged::Bool
    message::String
end

# Float32 数据使用 Float32 计算（节省内存，GPU 兼容）
xs32 = Float32[1.0, 2.0, 3.0, 4.0]
ys32 = Float32[1.0, 4.0, 9.0, 16.0]
interp32 = LinearInterp(xs32, ys32)  # LinearInterp{Float32}

# Float64 数据使用 Float64 计算
xs64 = Float64[1.0, 2.0, 3.0, 4.0]
ys64 = Float64[1.0, 4.0, 9.0, 16.0]
interp64 = LinearInterp(xs64, ys64)  # LinearInterp{Float64}

# BigFloat 高精度计算
xs_big = BigFloat[1, 2, 3, 4]
ys_big = BigFloat[1, 4, 9, 16]
interp_big = LinearInterp(xs_big, ys_big)  # LinearInterp{BigFloat}

# src/integration.jl

"""
    quad(f, a, b; n=1000) -> QuadratureResult

使用复合梯形法则计算 `f` 在 `[a, b]` 上的积分。

# 示例
```julia
result = quad(sin, 0, π)
result.value  # ≈ 2.0

# 高斯积分节点和权重（编译期常量）
const GAUSS_POINTS = Dict{Int, Vector{Float64}}(
    2 => [-0.5773502691896257, 0.5773502691896257],
    3 => [-0.7745966692414834, 0.0, 0.7745966692414834],
    4 => [-0.8611363115940526, -0.3399810435848563,
           0.3399810435848563,  0.8611363115940526],
    5 => [-0.9061798459386640, -0.5384693101056831, 0.0,
           0.5384693101056831,  0.9061798459386640],
)

const GAUSS_WEIGHTS = Dict{Int, Vector{Float64}}(
    2 => [1.0, 1.0],
    3 => [0.5555555555555556, 0.8888888888888888, 0.5555555555555556],
    4 => [0.3478548451374538, 0.6521451548625461,
          0.6521451548625461, 0.3478548451374538],
    5 => [0.2369268850561891, 0.4786286704993665, 0.5688888888888889,
          0.4786286704993665, 0.2369268850561891],
)

"""
    gauss_legendre(f, a, b; order=5)

使用 Gauss-Legendre 求积计算积分。对多项式精确。
"""
function gauss_legendre(f::F, a::T, b::T; order::Int=5) where {F, T<:Real}
    haskey(GAUSS_POINTS, order) || throw(ArgumentError("order 必须为 2-5"))
    
    points = GAUSS_POINTS[order]
    weights = GAUSS_WEIGHTS[order]
    
    # 映射到 [a, b]
    mid = (a + b) / 2
    half = (b - a) / 2
    
    s = zero(T)
    for (p, w) in zip(points, weights)
        s += w * f(mid + half * p)
    end
    
    return QuadratureResult{T}(s * half, zero(T), order, true)
end

# src/differentiation.jl

"""
    central_diff(f, x; h=1e-8)

使用中心差分计算导数。
"""
function central_diff(f::F, x::T; h::Real=1e-8) where {F, T<:Real}
    return (f(x + h) - f(x - h)) / (2h)
end

"""
    gradient(f, x; h=1e-8)

计算标量函数的梯度向量。
"""
function gradient(f::F, x::Vector{T}; h::Real=1e-8) where {F, T<:Real}
    n = length(x)
    grad = zeros(T, n)
    fx = f(x)
    
    for i in 1:n
        x_plus = copy(x)
        x_minus = copy(x)
        x_plus[i] += h
        x_minus[i] -= h
        grad[i] = (f(x_plus) - f(x_minus)) / (2h)
    end
    
    return grad
end

"""
    jacobian(f, x; h=1e-8)

计算向量函数的 Jacobian 矩阵。
"""
function jacobian(f::F, x::Vector{T}; h::Real=1e-8) where {F, T<:Real}
    n = length(x)
    fx = f(x)
    m = length(fx)
    J = zeros(T, m, n)
    
    for i in 1:n
        x_plus = copy(x)
        x_minus = copy(x)
        x_plus[i] += h
        x_minus[i] -= h
        J[:, i] = (f(x_plus) - f(x_minus)) / (2h)
    end
    
    return J
end

# src/interpolation.jl

"""
    (interp::LinearInterp)(x)

线性插值求值。
"""
function (interp::LinearInterp{T})(x::Real) where T
    x = T(x)
    
    # 边界检查
    (interp.xs[1] <= x <= interp.xs[end]) ||
        throw(DomainError(x, "超出插值范围"))
    
    # 二分查找
    lo, hi = 1, interp.n
    while hi - lo > 1
        mid = (lo + hi) ÷ 2
        if interp.xs[mid] <= x
            lo = mid
        else
            hi = mid
        end
    end
    
    # 线性插值
    t = (x - interp.xs[lo]) / (interp.xs[hi] - interp.xs[lo])
    return interp.ys[lo] + t * (interp.ys[hi] - interp.ys[lo])
end

"""
    CubicSpline(xs, ys)

构造三次样条插值器（自然边界条件）。
"""
function CubicSpline(xs::Vector{T}, ys::Vector{T}) where T
    n = length(xs)
    n >= 3 || throw(ArgumentError("至少需要 3 个数据点"))
    
    # 计算二阶导数（三对角系统）
    h = diff(xs)
    alpha = zeros(T, n)
    for i in 2:n-1
        alpha[i] = 3/h[i] * (ys[i+1] - ys[i]) - 3/h[i-1] * (ys[i] - ys[i-1])
    end
    
    # Thomas 算法求解三对角系统
    l = ones(T, n)
    mu = zeros(T, n)
    z = zeros(T, n)
    
    for i in 2:n-1
        l[i] = 2(xs[i+1] - xs[i-1]) - h[i-1] * mu[i-1]
        mu[i] = h[i] / l[i]
        z[i] = (alpha[i] - h[i-1] * z[i-1]) / l[i]
    end
    
    # 反向求解系数
    c = zeros(T, n)
    b = zeros(T, n-1)
    d = zeros(T, n-1)
    
    for j in n-1:-1:1
        c[j] = z[j] - mu[j] * c[j+1]
        b[j] = (ys[j+1] - ys[j]) / h[j] - h[j] * (c[j+1] + 2c[j]) / 3
        d[j] = (c[j+1] - c[j]) / (3h[j])
    end
    
    # 存储系数 [a, b, c, d]
    coeffs = zeros(T, n-1, 4)
    coeffs[:, 1] = ys[1:end-1]  # a
    coeffs[:, 2] = b             # b
    coeffs[:, 3] = c[1:end-1]   # c
    coeffs[:, 4] = d             # d
    
    return CubicSpline{T}(xs, ys, coeffs, n-1)
end

function (cs::CubicSpline{T})(x::Real) where T
    x = T(x)
    
    # 二分查找区间
    lo, hi = 1, cs.n
    while hi - lo > 1
        mid = (lo + hi) ÷ 2
        if cs.xs[mid] <= x
            lo = mid
        else
            hi = mid
        end
    end
    
    # 计算样条值
    dx = x - cs.xs[lo]
    c = cs.coeffs[lo, :]
    return c[1] + c[2]*dx + c[3]*dx^2 + c[4]*dx^3
end

# src/optimization.jl

"""
    gradient_descent(f, x0; lr=0.01, max_iter=1000, tol=1e-6)

梯度下降优化。
"""
function gradient_descent(f::F, x0::Vector{T};
                          lr::Real=0.01, max_iter::Int=1000,
                          tol::Real=1e-6) where {F, T<:Real}
    x = copy(x0)
    fx = f(x)
    
    for iter in 1:max_iter
        g = gradient(f, x; h=1e-8)
        x_new = x .- lr .* g
        fx_new = f(x_new)
        
        if norm(x_new - x) < tol
            return OptimizationResult{T}(x_new, fx_new, iter, true, "收敛")
        end
        
        x = x_new
        fx = fx_new
    end
    
    return OptimizationResult{T}(x, fx, max_iter, false, "达到最大迭代次数")
end

"""
    newton_method(f, x0; max_iter=100, tol=1e-8)

牛顿法优化（使用数值 Hessian）。
"""
function newton_method(f::F, x0::Vector{T};
                       max_iter::Int=100,
                       tol::Real=1e-8) where {F, T<:Real}
    x = copy(x0)
    n = length(x)
    
    for iter in 1:max_iter
        g = gradient(f, x; h=1e-6)
        
        if norm(g) < tol
            return OptimizationResult{T}(x, f(x), iter, true, "梯度接近零")
        end
        
        # 数值 Hessian
        H = zeros(T, n, n)
        h = 1e-5
        for i in 1:n
            for j in i:n
                xpp = copy(x); xpp[i] += h; xpp[j] += h
                xpm = copy(x); xpm[i] += h; xpm[j] -= h
                xmp = copy(x); xmp[i] -= h; xmp[j] += h
                xmm = copy(x); xmm[i] -= h; xmm[j] -= h
                H[i,j] = (f(xpp) - f(xpm) - f(xmp) + f(xmm)) / (4h^2)
                H[j,i] = H[i,j]
            end
        end
        
        # 加正则化确保正定
        H += 1e-6 * I
        
        x = x - H \ g
    end
    
    return OptimizationResult{T}(x, f(x), max_iter, false, "达到最大迭代次数")
end

# src/utils.jl

"""
    linspace(a, b, n)

生成等间距向量，类似 MATLAB 的 linspace。
"""
linspace(a, b, n) = range(a, b; length=n)

"""
    numdiff_table(f, x; steps)

打印数值微分的收敛表（Richardson 外推）。
"""
function numdiff_table(f, x; steps=[1.0/2^k for k in 0:10])
    exact = central_diff(f, x; h=1e-12)
    
    println("h\t\t估计值\t\t误差\t\t收敛比")
    println("-"^60)
    
    prev_est = nothing
    for h in steps
        est = central_diff(f, x; h=h)
        err = abs(est - exact)
        
        ratio = prev_est !== nothing ? abs(prev_est - exact) / max(err, 1e-16) : NaN
        
        @printf "%.2e\t%.8f\t%.2e\t%.2f\n" h est err ratio
        prev_est = est
    end
end

# test/test_integration.jl
using Test
using SciKit

@testset "数值积分" begin
    @testset "quad — 梯形法则" begin
        result = quad(sin, 0, π)
        @test result.value ≈ 2.0 atol=1e-6
        @test result.converged
        
        result = quad(x -> x^2, 0, 1)
        @test result.value ≈ 1/3 atol=1e-6
        
        # 零长度区间
        result = quad(sin, 1, 1)
        @test result.value ≈ 0.0
    end
    
    @testset "simpson — Simpson 法则" begin
        result = simpson(sin, 0, π)
        @test result.value ≈ 2.0 atol=1e-10
        
        # 多项式精确
        result = simpson(x -> x^3, 0, 1)
        @test result.value ≈ 0.25 atol=1e-12
    end
    
    @testset "gauss_legendre — 高斯积分" begin
        for order in 2:5
            result = gauss_legendre(x -> x^(2*order-2), -1, 1; order=order)
            @test result.value ≈ 2/(2*order-1) atol=1e-10
        end
    end
    
    @testset "参数验证" begin
        @test_throws ArgumentError quad(sin, 0, 1; n=0)
    end
end

# test/test_properties.jl
using Test
using SciKit

@testset "属性测试" begin
    @testset "线性性质" begin
        f(x) = 3x + 2
        result = quad(f, 0, 5)
        @test result.value ≈ 42.5  # ∫₀⁵ (3x+2) dx = 3·25/2 + 10 = 47.5
    end
    
    @testset "插值边界条件" begin
        xs = [0.0, 1.0, 2.0, 3.0]
        ys = [0.0, 1.0, 4.0, 9.0]
        interp = LinearInterp(xs, ys)
        
        # 插值器应通过所有数据点
        for (x, y) in zip(xs, ys)
            @test interp(x) ≈ y
        end
    end
    
    @testset "微分近似精度" begin
        # 中心差分误差应为 O(h²)
        f(x) = sin(x)
        exact = cos(1.0)
        
        err1 = abs(central_diff(f, 1.0; h=1e-2) - exact)
        err2 = abs(central_diff(f, 1.0; h=1e-4) - exact)
        
        # err1/err2 ≈ (1e-2/1e-4)² = 10000
        @test err1 / err2 > 100
    end
end

# test/test_integration_e2e.jl
@testset "端到端测试" begin
    # 综合使用积分+插值+微分
    xs = range(0, 2π, length=100)
    ys = sin.(xs)
    
    # 插值
    interp = LinearInterp(collect(xs), collect(ys))
    
    # 积分插值后的函数
    result = quad(interp, 0, 2π; n=10000)
    @test result.value ≈ 0.0 atol=1e-3  # sin 的完整周期积分为 0
end

# docs/make.jl
using Documenter
using SciKit

makedocs(;
    modules=[SciKit],
    sitename="SciKit.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://user.github.io/SciKit.jl",
        assets=["assets/custom.css"],
    ),
    pages=[
        "首页" => "index.md",
        "快速入门" => "getting-started.md",
        "教程" => [
            "数值积分" => "tutorials/integration.md",
            "数值微分" => "tutorials/differentiation.md",
            "插值" => "tutorials/interpolation.md",
        ],
        "API 参考" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/user/SciKit.jl",
    devbranch="main",
)

"""
    quad(f, a::Real, b::Real; n::Int=1000) -> QuadratureResult

使用复合梯形法则计算函数 `f` 在区间 `[a, b]` 上的定积分。

# 参数
- `f`: 被积函数，接受单个实数参数
- `a`: 积分下限
- `b`: 积分上限
- `n`: 子区间数量（默认 1000）

# 返回值
返回 `QuadratureResult`，包含：
- `value`: 积分近似值
- `error`: 误差估计
- `n_evals`: 函数求值次数
- `converged`: 是否收敛

# 示例
```julia
# 计算 ∫₀π sin(x) dx = 2
result = quad(sin, 0, π)
result.value  # ≈ 2.0

# 自定义函数
f(x) = exp(-x^2)
result = quad(f, -5, 5)  # 近似 √π

# .github/workflows/CI.yml
name: CI

on:
  push:
    branches: [main]
  pull_request:
    branches: [main]

jobs:
  test:
    name: Julia ${{ matrix.version }} - ${{ matrix.os }}
    runs-on: ${{ matrix.os }}
    strategy:
      matrix:
        version: ['1.10', '1.11']
        os: [ubuntu-latest, windows-latest, macOS-latest]
    steps:
      - uses: actions/checkout@v4
      - uses: julia-actions/setup-julia@v2
        with:
          version: ${{ matrix.version }}
      - uses: julia-actions/cache@v2
      - uses: julia-actions/julia-buildpkg@v1
      - uses: julia-actions/julia-runtest@v1
      - uses: julia-actions/julia-processcoverage@v1
      - uses: codecov/codecov-action@v4
        with:
          files: lcov.info

  docs:
    name: Documentation
    runs-on: ubuntu-latest
    steps:
      - uses: actions/checkout@v4
      - uses: julia-actions/setup-julia@v2
      - uses: julia-actions/cache@v2
      - run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
      - run: julia --project=docs/ docs/make.jl
        env:
          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

# 1. 确保测试通过
julia --project -e 'using Pkg; Pkg.test()'

# 2. 更新版本号
# Project.toml: version = "0.1.0"

# 3. 提交并推送
git add -A
git commit -m "SciKit.jl v0.1.0"
git tag v0.1.0
git push --tags

# 4. 注册
# 在 GitHub issue 中 @JuliaRegistrator register