#include "fiko/sketch/solver.hpp"

#include <algorithm>
#include <cmath>
#include <utility>

namespace fiko::eskiz {

namespace {

// ---- Yardimcilar: geometri kimliginden Geometri* + nokta secimi ----

const Geometri* geoBul(const std::vector<Geometri>& geo, u32 kimlik) {
  for (const auto& g : geo)
    if (g.kimlik == kimlik) return &g;
  return nullptr;
}

// Bir geometri + "hangi nokta" (0=a,1=b,2=merkez) icin Nokta* dondurur.
const Nokta* noktaSec(const Geometri& g, int hangi) {
  switch (hangi) {
    case 0: return g.a.gecerli() ? &g.a : (g.merkez.gecerli() ? &g.merkez : nullptr);
    case 1: return g.b.gecerli() ? &g.b : nullptr;
    case 2: return g.merkez.gecerli() ? &g.merkez : nullptr;
    default: return nullptr;
  }
}

struct N2 { Skalar x, y; };

N2 nk(const std::vector<Skalar>& p, const Nokta* n) {
  if (!n || !n->gecerli()) return {0.0, 0.0};
  return {p[(std::size_t)n->ix], p[(std::size_t)n->iy]};
}

// Cizgi yon vektoru (b - a).
struct Cizgi { N2 a, b; };

Cizgi cizgiAl(const std::vector<Skalar>& p, const Geometri& g) {
  return {nk(p, &g.a), nk(p, &g.b)};
}

// Bir geometrinin TUM param indeksleri (SPARSITY: kisit → ilgili Jacobian sutunlari;
// bir kisitin residual'i yalniz kendi geo'larinin parametrelerine baglidir).
void geoParamlari(const Geometri& g, std::vector<int>& out) {
  const auto n = [&](const Nokta& x) {
    if (x.ix >= 0) out.push_back(x.ix);
    if (x.iy >= 0) out.push_back(x.iy);
  };
  n(g.a); n(g.b); n(g.merkez);
  for (int i : {g.ir, g.ir2, g.iaci, g.ib1, g.ib2})
    if (i >= 0) out.push_back(i);
  for (const Nokta& e : g.ekNoktalar) n(e);
}

}  // namespace

int Solver::satirSayisi(const Kisit& k) {
  switch (k.tip) {
    case KisitTip::Cakisik:
    case KisitTip::Midpoint:
      return 2;  // dx, dy
    case KisitTip::Simetrik:
      return 2;  // orta-eksende + segment-eksene-dik (TAM simetri)
    case KisitTip::Sabit:
      return 0;  // residual uretmez; app param'i sabitParam ile dondurur
    default:
      return 1;
  }
}

namespace {
// TEK kisitin residual satir(lari)ni r'ye APPEND eder — artiklar + sparse Jacobian
// ORTAK cekirdegi (sparsity icin sart: perturbasyonda yalniz ilgili kisit yeniden
// hesaplanir, tum eskiz degil). Satir sayisi Solver::satirSayisi ile SOZLESMELI.
void kisitArtik(const std::vector<Skalar>& p, const std::vector<Geometri>& geo,
                const Kisit& k, std::vector<Skalar>& r) {
  // Ilgili geometrileri cek (guvenli — eksikse residual 0 birak).
  const Geometri* g0 = k.geo.size() > 0 ? geoBul(geo, k.geo[0]) : nullptr;
  const Geometri* g1 = k.geo.size() > 1 ? geoBul(geo, k.geo[1]) : nullptr;
  const int h0 = k.nokta.size() > 0 ? k.nokta[0] : 0;
  const int h1 = k.nokta.size() > 1 ? k.nokta[1] : 0;

  switch (k.tip) {
      case KisitTip::Cakisik: {
        // 2 nokta cakisik: f = (x2-x1, y2-y1)
        N2 p0 = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
        N2 p1 = nk(p, g1 ? noktaSec(*g1, h1) : nullptr);
        r.push_back(p1.x - p0.x);
        r.push_back(p1.y - p0.y);
        break;
      }
      case KisitTip::Sabit:
        // Residual uretmez (satirSayisi=0). App ilgili param indekslerini
        // coz() sabitParam kumesine koyarak dondurur.
        break;
      case KisitTip::Midpoint: {
        // Nokta (g0,h0) cizginin (g1) ORTA noktasinda: f = pt - (a+b)/2.
        N2 pt = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
        if (g1 && g1->tip == GeoTip::Cizgi) {
          Cizgi c = cizgiAl(p, *g1);
          r.push_back(pt.x - 0.5 * (c.a.x + c.b.x));
          r.push_back(pt.y - 0.5 * (c.a.y + c.b.y));
        } else {
          r.push_back(0.0);
          r.push_back(0.0);
        }
        break;
      }
      case KisitTip::NoktaUzerinde: {
        // g0 nokta (nokta h0), g1 cizgi veya daire.
        N2 pt = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
        if (g1 && g1->tip == GeoTip::Cizgi) {
          Cizgi c = cizgiAl(p, *g1);
          Skalar dx = c.b.x - c.a.x, dy = c.b.y - c.a.y;
          Skalar len = std::sqrt(dx * dx + dy * dy);
          // Isaretli dik mesafe: cross(d, pt-a)/|d|
          Skalar cross = dx * (pt.y - c.a.y) - dy * (pt.x - c.a.x);
          r.push_back(len > 1e-12 ? cross / len : cross);
        } else if (g1 && (g1->tip == GeoTip::Daire || g1->tip == GeoTip::Yay)) {
          N2 m = nk(p, &g1->merkez);
          Skalar rad = p[(std::size_t)g1->ir];
          Skalar dx = pt.x - m.x, dy = pt.y - m.y;
          r.push_back(std::sqrt(dx * dx + dy * dy) - rad);
        } else {
          r.push_back(0.0);
        }
        break;
      }
      case KisitTip::Yatay: {
        // Cizgi ise a,b uclari; iki nokta ise g0(h0)/g1(h1). f = y2 - y1
        if (g0 && g0->tip == GeoTip::Cizgi && k.geo.size() == 1) {
          Cizgi c = cizgiAl(p, *g0);
          r.push_back(c.b.y - c.a.y);
        } else {
          N2 a = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
          N2 b = nk(p, g1 ? noktaSec(*g1, h1) : nullptr);
          r.push_back(b.y - a.y);
        }
        break;
      }
      case KisitTip::Dikey: {
        if (g0 && g0->tip == GeoTip::Cizgi && k.geo.size() == 1) {
          Cizgi c = cizgiAl(p, *g0);
          r.push_back(c.b.x - c.a.x);
        } else {
          N2 a = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
          N2 b = nk(p, g1 ? noktaSec(*g1, h1) : nullptr);
          r.push_back(b.x - a.x);
        }
        break;
      }
      case KisitTip::Paralel: {
        // 2 cizgi: f = sin(aci) = (d0 × d1) / (|d0||d1|) = 0
        // Boyutsuz olculur ki residual O(1); aksi halde cizgi uzunluguyla
        // olceklenir JᵀJ kotu-kosullu olur ve LM cok yavas yakinsar.
        Cizgi c0 = g0 ? cizgiAl(p, *g0) : Cizgi{};
        Cizgi c1 = g1 ? cizgiAl(p, *g1) : Cizgi{};
        Skalar d0x = c0.b.x - c0.a.x, d0y = c0.b.y - c0.a.y;
        Skalar d1x = c1.b.x - c1.a.x, d1y = c1.b.y - c1.a.y;
        Skalar len0 = std::sqrt(d0x * d0x + d0y * d0y);
        Skalar len1 = std::sqrt(d1x * d1x + d1y * d1y);
        Skalar olcek = len0 * len1;
        Skalar cross = d0x * d1y - d0y * d1x;
        r.push_back(olcek > 1e-12 ? cross / olcek : cross);
        break;
      }
      case KisitTip::Dik: {
        // 2 cizgi: f = cos(aci) = (d0 · d1) / (|d0||d1|) = 0
        // Boyutsuz olculur (yukaridaki Paralel notu gecerli).
        Cizgi c0 = g0 ? cizgiAl(p, *g0) : Cizgi{};
        Cizgi c1 = g1 ? cizgiAl(p, *g1) : Cizgi{};
        Skalar d0x = c0.b.x - c0.a.x, d0y = c0.b.y - c0.a.y;
        Skalar d1x = c1.b.x - c1.a.x, d1y = c1.b.y - c1.a.y;
        Skalar len0 = std::sqrt(d0x * d0x + d0y * d0y);
        Skalar len1 = std::sqrt(d1x * d1x + d1y * d1y);
        Skalar olcek = len0 * len1;
        Skalar dot = d0x * d1x + d0y * d1y;
        r.push_back(olcek > 1e-12 ? dot / olcek : dot);
        break;
      }
      case KisitTip::Tegel: {
        // cizgi-daire: |merkez-cizgi dik mesafe| - r = 0
        // daire-daire: |m1-m2| - (r1+r2) = 0  (dis teget)
        if (g0 && g1 && g0->tip == GeoTip::Cizgi &&
            (g1->tip == GeoTip::Daire || g1->tip == GeoTip::Yay)) {
          Cizgi c = cizgiAl(p, *g0);
          N2 m = nk(p, &g1->merkez);
          Skalar rad = p[(std::size_t)g1->ir];
          Skalar dx = c.b.x - c.a.x, dy = c.b.y - c.a.y;
          Skalar len = std::sqrt(dx * dx + dy * dy);
          Skalar cross = dx * (m.y - c.a.y) - dy * (m.x - c.a.x);
          Skalar dist = len > 1e-12 ? std::abs(cross / len) : std::abs(cross);
          r.push_back(dist - rad);
        } else if (g0 && g1 &&
                   (g0->tip == GeoTip::Daire || g0->tip == GeoTip::Yay) &&
                   (g1->tip == GeoTip::Daire || g1->tip == GeoTip::Yay)) {
          N2 m0 = nk(p, &g0->merkez), m1 = nk(p, &g1->merkez);
          Skalar r0 = p[(std::size_t)g0->ir], r1 = p[(std::size_t)g1->ir];
          Skalar dx = m1.x - m0.x, dy = m1.y - m0.y;
          r.push_back(std::sqrt(dx * dx + dy * dy) - (r0 + r1));
        } else {
          r.push_back(0.0);
        }
        break;
      }
      case KisitTip::Esit: {
        // 2 cizgi: esit uzunluk. 2 daire/yay: esit yaricap.
        if (g0 && g1 && g0->tip == GeoTip::Cizgi && g1->tip == GeoTip::Cizgi) {
          Cizgi c0 = cizgiAl(p, *g0), c1 = cizgiAl(p, *g1);
          Skalar l0 = std::sqrt((c0.b.x - c0.a.x) * (c0.b.x - c0.a.x) +
                                (c0.b.y - c0.a.y) * (c0.b.y - c0.a.y));
          Skalar l1 = std::sqrt((c1.b.x - c1.a.x) * (c1.b.x - c1.a.x) +
                                (c1.b.y - c1.a.y) * (c1.b.y - c1.a.y));
          r.push_back(l1 - l0);
        } else if (g0 && g1 && g0->ir >= 0 && g1->ir >= 0) {
          r.push_back(p[(std::size_t)g1->ir] - p[(std::size_t)g0->ir]);
        } else {
          r.push_back(0.0);
        }
        break;
      }
      case KisitTip::Simetrik: {
        // g0,g1 iki nokta; g2 simetri cizgisi. TAM simetri = 2 kosul:
        //   (1) orta nokta eksen dogrusunda (dik mesafe = 0)
        //   (2) a→b segmenti eksene DIK (dot = 0)
        // Eski kod yalniz (1)'i uyguluyordu ("yarim") → noktalar eksen boyunca
        // kayabiliyordu; simetri acisi serbest kaliyordu. 07-22 tamamlandi.
        const Geometri* g2 = k.geo.size() > 2 ? geoBul(geo, k.geo[2]) : nullptr;
        N2 a = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
        N2 b = nk(p, g1 ? noktaSec(*g1, h1) : nullptr);
        if (g2 && g2->tip == GeoTip::Cizgi) {
          Cizgi c = cizgiAl(p, *g2);
          N2 orta{(a.x + b.x) * 0.5, (a.y + b.y) * 0.5};
          Skalar dx = c.b.x - c.a.x, dy = c.b.y - c.a.y;
          Skalar len = std::sqrt(dx * dx + dy * dy);
          Skalar cross = dx * (orta.y - c.a.y) - dy * (orta.x - c.a.x);
          Skalar dik = dx * (b.x - a.x) + dy * (b.y - a.y);
          r.push_back(len > 1e-12 ? cross / len : cross);
          r.push_back(len > 1e-12 ? dik / len : dik);
        } else {
          r.push_back(0.0);
          r.push_back(0.0);
        }
        break;
      }
      case KisitTip::Mesafe: {
        // Cizgi (tek geo): uzunluk - deger.  Iki nokta: mesafe - deger.
        if (g0 && g0->tip == GeoTip::Cizgi && k.geo.size() == 1) {
          Cizgi c = cizgiAl(p, *g0);
          Skalar dx = c.b.x - c.a.x, dy = c.b.y - c.a.y;
          r.push_back(std::sqrt(dx * dx + dy * dy) - k.deger);
        } else {
          N2 a = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
          N2 b = nk(p, g1 ? noktaSec(*g1, h1) : nullptr);
          Skalar dx = b.x - a.x, dy = b.y - a.y;
          r.push_back(std::sqrt(dx * dx + dy * dy) - k.deger);
        }
        break;
      }
      case KisitTip::YatayMesafe: {
        if (g0 && g0->tip == GeoTip::Cizgi && k.geo.size() == 1) {
          Cizgi c = cizgiAl(p, *g0);
          r.push_back(c.b.x - c.a.x - k.deger);
        } else {
          N2 a = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
          N2 b = nk(p, g1 ? noktaSec(*g1, h1) : nullptr);
          r.push_back(b.x - a.x - k.deger);
        }
        break;
      }
      case KisitTip::DikeyMesafe: {
        if (g0 && g0->tip == GeoTip::Cizgi && k.geo.size() == 1) {
          Cizgi c = cizgiAl(p, *g0);
          r.push_back(c.b.y - c.a.y - k.deger);
        } else {
          N2 a = nk(p, g0 ? noktaSec(*g0, h0) : nullptr);
          N2 b = nk(p, g1 ? noktaSec(*g1, h1) : nullptr);
          r.push_back(b.y - a.y - k.deger);
        }
        break;
      }
      case KisitTip::Yaricap: {
        if (g0 && g0->ir >= 0)
          r.push_back(p[(std::size_t)g0->ir] - k.deger);
        else
          r.push_back(0.0);
        break;
      }
      case KisitTip::Cap: {
        if (g0 && g0->ir >= 0)
          r.push_back(2.0 * p[(std::size_t)g0->ir] - k.deger);
        else
          r.push_back(0.0);
        break;
      }
      case KisitTip::Aci: {
        // 2 cizgi arasi aci = deger (radyan). f = atan2(cross,dot) - deger
        Cizgi c0 = g0 ? cizgiAl(p, *g0) : Cizgi{};
        Cizgi c1 = g1 ? cizgiAl(p, *g1) : Cizgi{};
        Skalar d0x = c0.b.x - c0.a.x, d0y = c0.b.y - c0.a.y;
        Skalar d1x = c1.b.x - c1.a.x, d1y = c1.b.y - c1.a.y;
        Skalar cross = d0x * d1y - d0y * d1x;
        Skalar dot = d0x * d1x + d0y * d1y;
        r.push_back(std::atan2(cross, dot) - k.deger);
        break;
      }
      case KisitTip::Semimajor: {
        // Elips buyuk yaricap (a) = deger. ir buyuk yaricap param indeksi.
        if (g0 && g0->ir >= 0)
          r.push_back(p[(std::size_t)g0->ir] - k.deger);
        else
          r.push_back(0.0);
        break;
      }
      case KisitTip::Semiminor: {
        // Elips kucuk yaricap (b) = deger. ir2 kucuk yaricap param indeksi.
        if (g0 && g0->ir2 >= 0)
          r.push_back(p[(std::size_t)g0->ir2] - k.deger);
        else
          r.push_back(0.0);
        break;
      }
    }
}
}  // namespace

// residual: her kisit satir(lari)ni r'ye append eder. Sira coz() ile ayni.
void Solver::artiklar(const std::vector<Skalar>& p,
                      const std::vector<Geometri>& geo,
                      const std::vector<Kisit>& kisit, std::vector<Skalar>& r) {
  r.clear();
  for (const auto& k : kisit) {
    if (!k.suru) continue;  // reference kisit cozumu surmez
    kisitArtik(p, geo, k, r);
  }
}

namespace {

// Dense NxN lineer cozum: Gauss eliminasyonu + kismi pivot. A|b -> x.
// Tekil (kucuk pivot) durumunda false doner (LM disardan λ artirir).
bool gaussCoz(std::vector<std::vector<Skalar>>& A, std::vector<Skalar>& b,
              std::vector<Skalar>& x, Skalar eps = 1e-14) {
  const int n = (int)b.size();
  for (int i = 0; i < n; ++i) {
    // pivot sec
    int pr = i;
    Skalar en = std::abs(A[(std::size_t)i][(std::size_t)i]);
    for (int r2 = i + 1; r2 < n; ++r2) {
      Skalar v = std::abs(A[(std::size_t)r2][(std::size_t)i]);
      if (v > en) { en = v; pr = r2; }
    }
    if (en < eps) return false;
    if (pr != i) {
      std::swap(A[(std::size_t)i], A[(std::size_t)pr]);
      std::swap(b[(std::size_t)i], b[(std::size_t)pr]);
    }
    // elimine
    for (int r2 = i + 1; r2 < n; ++r2) {
      Skalar f = A[(std::size_t)r2][(std::size_t)i] / A[(std::size_t)i][(std::size_t)i];
      if (f == 0.0) continue;
      for (int c = i; c < n; ++c)
        A[(std::size_t)r2][(std::size_t)c] -= f * A[(std::size_t)i][(std::size_t)c];
      b[(std::size_t)r2] -= f * b[(std::size_t)i];
    }
  }
  x.assign((std::size_t)n, 0.0);
  for (int i = n - 1; i >= 0; --i) {
    Skalar s = b[(std::size_t)i];
    for (int c = i + 1; c < n; ++c)
      s -= A[(std::size_t)i][(std::size_t)c] * x[(std::size_t)c];
    x[(std::size_t)i] = s / A[(std::size_t)i][(std::size_t)i];
  }
  return true;
}

// Matris rank'i satir indirgeme (Gauss, kismi pivot) ile. Serbest parametre
// sutunlari sabitse (fixed) o sutunlar J'den cikarilmis olarak gelmelidir.
int matrisRank(std::vector<std::vector<Skalar>> M, int satir, int sutun,
               std::vector<int>* pivotsuz = nullptr, Skalar eps = 1e-9) {
  int rank = 0;
  for (int c = 0; c < sutun; ++c) {
    if (rank >= satir) {  // tüm satırlar pivotlandı → kalan sütunlar serbest
      if (pivotsuz) pivotsuz->push_back(c);
      continue;
    }
    // pivot bul
    int pr = -1;
    Skalar en = eps;
    for (int r2 = rank; r2 < satir; ++r2) {
      Skalar v = std::abs(M[(std::size_t)r2][(std::size_t)c]);
      if (v > en) { en = v; pr = r2; }
    }
    if (pr < 0) { if (pivotsuz) pivotsuz->push_back(c); continue; }
    std::swap(M[(std::size_t)rank], M[(std::size_t)pr]);
    Skalar piv = M[(std::size_t)rank][(std::size_t)c];
    for (int r2 = 0; r2 < satir; ++r2) {
      if (r2 == rank) continue;
      Skalar f = M[(std::size_t)r2][(std::size_t)c] / piv;
      if (f == 0.0) continue;
      for (int cc = c; cc < sutun; ++cc)
        M[(std::size_t)r2][(std::size_t)cc] -= f * M[(std::size_t)rank][(std::size_t)cc];
    }
    ++rank;
  }
  return rank;
}

Skalar normSqr(const std::vector<Skalar>& v) {
  Skalar s = 0;
  for (Skalar x : v) s += x * x;
  return s;
}

}  // namespace

Solver::Sonuc Solver::coz(std::vector<Skalar>& p,
                          const std::vector<Geometri>& geo,
                          const std::vector<Kisit>& kisit,
                          const std::set<int>* sabitParam) {
  Sonuc sonuc;

  // Serbest parametre indeksleri (sabit olmayanlar) → azalt-map.
  const int np = (int)p.size();
  std::vector<int> serbest;  // serbest sutun -> global param indeks
  for (int i = 0; i < np; ++i) {
    if (sabitParam && sabitParam->count(i)) continue;
    serbest.push_back(i);
  }
  const int nf = (int)serbest.size();

  // Suru (driving) kisit satir sayisi.
  int m = 0;
  for (const auto& k : kisit)
    if (k.suru) m += satirSayisi(k);

  // Kisit yoksa: sadece dof hesapla, cozme. Tum serbest paramlar kisitsiz.
  if (m == 0) {
    sonuc.dof = nf;
    sonuc.durum = nf == 0 ? Durum::IyiKisitli : Durum::AzKisitli;
    sonuc.kalanNorm = 0;
    sonuc.serbestParamlar = serbest;
    return sonuc;
  }

  std::vector<Skalar> r;
  artiklar(p, geo, kisit, r);
  Skalar hataSqr = normSqr(r);

  // Jacobian buffer (merkezi sonlu-fark).
  std::vector<std::vector<Skalar>> J(
      (std::size_t)m, std::vector<Skalar>((std::size_t)nf, 0.0));

  // ---- SPARSITY (ölçüldü: 60 çizgi/119 kısıt sürükle 183 ms — tam-yoğun fark her
  // satır için TÜM nf paramı pertürbe ediyordu). Bir kısıtın residual'ı yalnız kendi
  // geo'larının parametrelerine bağlı → J yalnız (kısıt-satırı, ilgili-param)
  // hücrelerinde sonlu-farkla doldurulur; kalanı SIFIR (matematiksel olarak da sıfır
  // — davranış birebir). Pertürbasyon tek kısıtı yeniden hesaplar (kisitArtik).
  std::vector<int> sutunMap((std::size_t)np, -1);
  for (int c = 0; c < nf; ++c) sutunMap[(std::size_t)serbest[(std::size_t)c]] = c;
  struct KBlok {
    const Kisit* k;
    int ro, satir;
    std::vector<std::pair<int, int>> ciftler;  // (param indeksi, serbest sütun)
  };
  std::vector<KBlok> bloklar;
  {
    int ro = 0;
    std::vector<int> pr;
    for (const auto& k : kisit) {
      if (!k.suru) continue;
      KBlok b;
      b.k = &k;
      b.ro = ro;
      b.satir = satirSayisi(k);
      ro += b.satir;
      if (b.satir == 0) continue;  // Sabit: satır üretmez
      pr.clear();
      for (u32 gk : k.geo)
        if (const Geometri* g = geoBul(geo, gk)) geoParamlari(*g, pr);
      std::sort(pr.begin(), pr.end());
      pr.erase(std::unique(pr.begin(), pr.end()), pr.end());
      for (int gi : pr) {
        const int c = (gi >= 0 && gi < np) ? sutunMap[(std::size_t)gi] : -1;
        if (c >= 0) b.ciftler.emplace_back(gi, c);
      }
      bloklar.push_back(std::move(b));
    }
  }
  std::vector<Skalar> bp, bm;
  // Satir basina nonzero listesi (sutun, deger) — A=JᵀJ ve Jᵀr birikimi bunlardan
  // (faz 2: yogun A hesabi O(nf²·m)≈110M carpim/iter baskin cikmisti; nonzero'dan
  // birikim satir basina ~8² carpima duser). Dense J YALNIZ rank icin de dolar (ucuz).
  std::vector<std::vector<std::pair<int, Skalar>>> Jsp((std::size_t)m);
  const auto jacobianKur = [&]() {
    for (auto& satirVec : J) std::fill(satirVec.begin(), satirVec.end(), 0.0);
    for (auto& sv : Jsp) sv.clear();
    for (const KBlok& b : bloklar) {
      for (const auto& [gi, c] : b.ciftler) {
        const Skalar eski = p[(std::size_t)gi];
        const Skalar h = fdAdim * (std::abs(eski) > 1.0 ? std::abs(eski) : 1.0);
        p[(std::size_t)gi] = eski + h;
        bp.clear();
        kisitArtik(p, geo, *b.k, bp);
        p[(std::size_t)gi] = eski - h;
        bm.clear();
        kisitArtik(p, geo, *b.k, bm);
        p[(std::size_t)gi] = eski;
        const Skalar inv = 1.0 / (2.0 * h);
        const int n2 = std::min(b.satir, (int)std::min(bp.size(), bm.size()));
        for (int j = 0; j < n2; ++j) {
          const Skalar v = (bp[(std::size_t)j] - bm[(std::size_t)j]) * inv;
          J[(std::size_t)(b.ro + j)][(std::size_t)c] = v;
          if (v != 0.0) Jsp[(std::size_t)(b.ro + j)].emplace_back(c, v);
        }
      }
    }
  };

  Skalar lambda = lambdaBas;
  int iter = 0;
  for (; iter < maxIter; ++iter) {
    if (std::sqrt(hataSqr) < tol) break;

    // --- Jacobian: SPARSE merkezi fark (yalniz ilgili kisit×param hücreleri) ---
    jacobianKur();

    // --- Normal denklemler: (JᵀJ + λ diag(JᵀJ)) Δ = −Jᵀr ---
    // NONZERO birikimi (faz 2): yogun O(nf²·m) yerine satir basina nonzero² —
    // matematiksel olarak AYNI A/g (sifir katkilar atlanir, davranis birebir).
    std::vector<std::vector<Skalar>> A(
        (std::size_t)nf, std::vector<Skalar>((std::size_t)nf, 0.0));
    std::vector<Skalar> g((std::size_t)nf, 0.0);
    for (int rr = 0; rr < m; ++rr) {
      const Skalar rv = r[(std::size_t)rr];
      for (const auto& [ci, vi] : Jsp[(std::size_t)rr]) {
        for (const auto& [cj, vj] : Jsp[(std::size_t)rr])
          A[(std::size_t)ci][(std::size_t)cj] += vi * vj;
        g[(std::size_t)ci] -= vi * rv;  // −Jᵀr
      }
    }

    // LM ic dongusu: λ ayarlayarak kabul edilen adim bul.
    bool adimKabul = false;
    for (int lmDeneme = 0; lmDeneme < 12 && !adimKabul; ++lmDeneme) {
      std::vector<std::vector<Skalar>> Ad = A;
      for (int i = 0; i < nf; ++i)
        Ad[(std::size_t)i][(std::size_t)i] += lambda * (A[(std::size_t)i][(std::size_t)i] + 1e-12);

      std::vector<Skalar> gd = g, delta;
      if (!gaussCoz(Ad, gd, delta)) {
        lambda *= 10.0;
        continue;
      }

      // Aday p.
      std::vector<Skalar> aday = p;
      for (int c = 0; c < nf; ++c)
        aday[(std::size_t)serbest[(std::size_t)c]] += delta[(std::size_t)c];

      std::vector<Skalar> rAday;
      artiklar(aday, geo, kisit, rAday);
      Skalar yeniSqr = normSqr(rAday);

      if (yeniSqr < hataSqr) {
        p = std::move(aday);
        r = std::move(rAday);
        hataSqr = yeniSqr;
        lambda = std::max(lambda * 0.3, 1e-12);
        adimKabul = true;
      } else {
        lambda *= 10.0;
      }
    }

    if (!adimKabul) break;  // iyilesme yok → dur (celiski veya minimum)
  }

  sonuc.iterasyon = iter;
  sonuc.kalanNorm = std::sqrt(hataSqr);

  // --- DOF + durum: son Jacobian'in rank'i ---
  // (Jacobian son iterasyondaki degerle; yakinsak noktaya yakin degerlendir.)
  jacobianKur();  // sparse — yoğun döngüyle aynı değerler (sıfırlar hesaplanmaz)
  std::vector<int> pivotsuz;
  const int rank = matrisRank(J, m, nf, &pivotsuz);
  sonuc.dof = nf - rank;
  sonuc.serbestParamlar.reserve(pivotsuz.size());
  for (int c : pivotsuz) sonuc.serbestParamlar.push_back(serbest[(std::size_t)c]);

  const bool yakinsadi = sonuc.kalanNorm < 1e-6;
  if (!yakinsadi) {
    // Yakinsamadi: fazla kisit varsa celiski, yoksa sadece yakinsamadi.
    sonuc.durum = (rank < m) ? Durum::Celiskili : Durum::Yakinsamadi;
  } else if (rank < m) {
    // Fazla kisit ama tutarli (redundant).
    sonuc.durum = Durum::AsiriKisitli;
  } else if (sonuc.dof == 0) {
    sonuc.durum = Durum::IyiKisitli;
  } else {
    sonuc.durum = Durum::AzKisitli;
  }
  return sonuc;
}

Solver::Sonuc Solver::surukle(std::vector<Skalar>& p,
                              const std::vector<Geometri>& geo,
                              const std::vector<Kisit>& kisit, u32 noktaKimlik,
                              Skalar hedefX, Skalar hedefY) {
  // Gecici yumusak kisit: noktaKimlik geometrisinin 'a' noktasini hedefe cek.
  const Geometri* g = geoBul(geo, noktaKimlik);
  if (!g) {
    return coz(p, geo, kisit);
  }
  const Nokta* n = noktaSec(*g, 0);
  if (!n || !n->gecerli()) {
    return coz(p, geo, kisit);
  }

  // Baslangic tahmini: noktayi dogrudan hedefe koy, sonra kisitlarla coz.
  // Gercek surukleme kisiti eklemeden — mevcut kisitlar noktayi geri ceker,
  // capa olarak hedefe yakin cozum bulunur. Daha guclu cekim icin gecici
  // YatayMesafe/DikeyMesafe kisitlari eklenebilir; burada dogrudan set + coz.
  std::vector<Kisit> genis = kisit;

  // Gecici bir "nokta noktada" hedefi: yeni bir sanal geo eklemek yerine
  // sabit hedef noktaya iki mesafe kisiti ekle (x ve y). Sanal capa noktasini
  // p'ye ekleyip sabit tut.
  std::vector<Skalar> pp = p;
  const int hx = (int)pp.size();
  pp.push_back(hedefX);
  const int hy = (int)pp.size();
  pp.push_back(hedefY);

  // Gecici geometri: sanal capa noktasi (kimlik cakismayan buyuk deger).
  std::vector<Geometri> genisGeo = geo;
  Geometri capa;
  capa.tip = GeoTip::Nokta;
  capa.a = Nokta{hx, hy};
  capa.kimlik = 0xFFFFFFFFu;
  genisGeo.push_back(capa);

  // Cakisik(surukle noktasi, capa) — capa parametreleri sabit.
  Kisit cek;
  cek.tip = KisitTip::Cakisik;
  cek.geo = {noktaKimlik, capa.kimlik};
  cek.nokta = {0, 0};
  cek.suru = true;
  cek.kimlik = 0xFFFFFFFEu;
  genis.push_back(cek);

  std::set<int> sabit{hx, hy};  // capa hedefte dondurulur
  Sonuc s = coz(pp, genisGeo, genis, &sabit);

  // Gercek parametreleri geri yaz (sanal capa haric).
  p.assign(pp.begin(), pp.begin() + (std::ptrdiff_t)p.size());
  return s;
}

}  // namespace fiko::eskiz