// OpenSTA, Static Timing Analyzer
// Copyright (c) 2026, Parallax Software, Inc.
// 
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
// 
// The origin of this software must not be misrepresented; you must not
// claim that you wrote the original software.
// 
// Altered source versions must be plainly marked as such, and must not be
// misrepresented as being the original software.
// 
// This notice may not be removed or altered from any source distribution.

#include "ReduceParasitics.hh"

#include <algorithm>
#include <map>
#include <set>

#include "Debug.hh"
#include "Error.hh"
#include "Liberty.hh"
#include "MinMax.hh"
#include "Network.hh"
#include "Parasitics.hh"
#include "Scene.hh"
#include "Sdc.hh"

namespace sta {

using ParasiticNodeValueMap = std::map<ParasiticNode *, double>;
using ResistorCurrentMap = std::map<ParasiticResistor *, double>;
using ParasiticResistorSet = std::set<ParasiticResistor *>;
using ParasiticNodeSet = std::set<ParasiticNode *>;

class ReduceToPi : public StaState
{
public:
  ReduceToPi(StaState *sta);
  void reduceToPi(const Parasitic *parasitic_network,
                  const Pin *drvr_pin,
                  ParasiticNode *drvr_node,
                  float coupling_cap_factor,
                  const RiseFall *rf,
                  const Scene *scene,
                  const MinMax *min_max,
                  float &c2,
                  float &rpi,
                  float &c1);
  bool pinCapsOneValue() { return pin_caps_one_value_; }

protected:
  void reducePiDfs(const Pin *drvr_pin,
                   ParasiticNode *node,
                   ParasiticResistor *from_res,
                   double src_resistance,
                   double &y1,
                   double &y2,
                   double &y3,
                   double &dwn_cap,
                   double &max_resistance);
  void visit(ParasiticNode *node);
  bool isVisited(ParasiticNode *node);
  void leave(ParasiticNode *node);
  void setDownstreamCap(ParasiticNode *node,
                        float cap);
  float downstreamCap(ParasiticNode *node);
  float pinCapacitance(ParasiticNode *node);
  bool isLoopResistor(ParasiticResistor *resistor);
  void markLoopResistor(ParasiticResistor *resistor);

  Parasitics *parasitics_;
  bool includes_pin_caps_;
  float coupling_cap_multiplier_ {1.0};
  const RiseFall *rf_ {nullptr};
  const Scene *scene_ {nullptr};
  const MinMax *min_max_ {nullptr};
  ParasiticNodeResistorMap resistor_map_;
  ParasiticNodeCapacitorMap capacitor_map_;

  ParasiticNodeSet visited_nodes_;
  ParasiticNodeValueMap node_values_;
  ParasiticResistorSet loop_resistors_;
  bool pin_caps_one_value_ {true};
};

ReduceToPi::ReduceToPi(StaState *sta) :
  StaState(sta)
{
}

// "Modeling the Driving-Point Characteristic of Resistive
// Interconnect for Accurate Delay Estimation", Peter O'Brien and
// Thomas Savarino, Proceedings of the 1989 Design Automation
// Conference.
void
ReduceToPi::reduceToPi(const Parasitic *parasitic_network,
                       const Pin *drvr_pin,
                       ParasiticNode *drvr_node,
                       float coupling_cap_factor,
                       const RiseFall *rf,
                       const Scene *scene,
                       const MinMax *min_max,
                       float &c2,
                       float &rpi,
                       float &c1)
{
  coupling_cap_multiplier_ = coupling_cap_factor;
  rf_ = rf;
  scene_ = scene;
  min_max_ = min_max;
  parasitics_ = scene_->parasitics(min_max);
  includes_pin_caps_ = parasitics_->includesPinCaps(parasitic_network),
  resistor_map_ = parasitics_->parasiticNodeResistorMap(parasitic_network);
  capacitor_map_ = parasitics_->parasiticNodeCapacitorMap(parasitic_network);

  double y1, y2, y3, dcap;
  double max_resistance = 0.0;
  reducePiDfs(drvr_pin, drvr_node, nullptr, 0.0,
              y1, y2, y3, dcap, max_resistance);

  if (y2 == 0.0 && y3 == 0.0) {
    // Capacitive load.
    c1 = y1;
    c2 = 0.0;
    rpi = 0.0;
  }
  else {
    c1 = y2 * y2 / y3;
    c2 = y1 - y2 * y2 / y3;
    rpi = -y3 * y3 / (y2 * y2 * y2);
  }
  debugPrint(debug_, "parasitic_reduce", 2,
             " Pi model c2={:.3g} rpi={:.3g} c1={:.3g} max_r={:.3g}",
             c2, rpi, c1, max_resistance);
}

// Find admittance moments.
void
ReduceToPi::reducePiDfs(const Pin *drvr_pin,
                        ParasiticNode *node,
                        ParasiticResistor *from_res,
                        double src_resistance,
                        double &y1,
                        double &y2,
                        double &y3,
                        double &dwn_cap,
                        double &max_resistance)
{
  double coupling_cap = 0.0;
  ParasiticCapacitorSeq &capacitors = capacitor_map_[node];
  for (ParasiticCapacitor *capacitor : capacitors)
    coupling_cap += parasitics_->value(capacitor);

  dwn_cap = parasitics_->nodeGndCap(node)
    + coupling_cap * coupling_cap_multiplier_
    + pinCapacitance(node);
  y1 = dwn_cap;
  y2 = y3 = 0.0;
  max_resistance = std::max(max_resistance, src_resistance);

  visit(node);
  ParasiticResistorSeq &resistors = resistor_map_[node];
  for (ParasiticResistor *resistor : resistors) {
    if (!isLoopResistor(resistor)) {
      ParasiticNode *onode = parasitics_->otherNode(resistor, node);
      // One commercial extractor creates resistors with identical from/to nodes.
      if (onode != node
          && resistor != from_res) {
        if (isVisited(onode)) {
          // Resistor loop.
          debugPrint(debug_, "parasitic_reduce", 2, " loop detected thru resistor {}",
                     parasitics_->id(resistor));
          markLoopResistor(resistor);
        }
        else {
          double r = parasitics_->value(resistor);
          double yd1, yd2, yd3, dcap;
          reducePiDfs(drvr_pin, onode, resistor, src_resistance + r,
                      yd1, yd2, yd3, dcap, max_resistance);
          // Rule 3.  Upstream traversal of a series resistor.
          // Rule 4.  Parallel admittances add.
          y1 += yd1;
          y2 += yd2 - r * yd1 * yd1;
          y3 += yd3 - 2 * r * yd1 * yd2 + r * r * yd1 * yd1 * yd1;
          dwn_cap += dcap;
        }
      }
    }
  }

  setDownstreamCap(node, dwn_cap);
  leave(node);
  debugPrint(debug_, "parasitic_reduce", 3,
             " node {} y1={:.3g} y2={:.3g} y3={:.3g} cap={:.3g}",
             parasitics_->name(node), y1, y2, y3, dwn_cap);
}

float
ReduceToPi::pinCapacitance(ParasiticNode *node)
{
  const Pin *pin = parasitics_->pin(node);
  float pin_cap = 0.0;
  if (pin) {
    Port *port = network_->port(pin);
    LibertyPort *lib_port = network_->libertyPort(port);
    const Sdc *sdc = scene_->sdc();
    if (lib_port) {
      if (!includes_pin_caps_) {
        pin_cap = sdc->pinCapacitance(pin, rf_, scene_, min_max_);
        pin_caps_one_value_ &= lib_port->capacitanceIsOneValue();
      }
    }
    else if (network_->isTopLevelPort(pin))
      pin_cap = sdc->portExtCap(port, rf_, min_max_);
  }
  return pin_cap;
}

void
ReduceToPi::visit(ParasiticNode *node)
{
  visited_nodes_.insert(node);
}

bool
ReduceToPi::isVisited(ParasiticNode *node)
{
  return visited_nodes_.contains(node);
}

void
ReduceToPi::leave(ParasiticNode *node)
{
  visited_nodes_.erase(node);
}

bool
ReduceToPi::isLoopResistor(ParasiticResistor *resistor)
{
  return loop_resistors_.contains(resistor);
}

void
ReduceToPi::markLoopResistor(ParasiticResistor *resistor)
{
  loop_resistors_.insert(resistor);
}

void
ReduceToPi::setDownstreamCap(ParasiticNode *node,
                             float cap)
{
  node_values_[node] = cap;
}

float
ReduceToPi::downstreamCap(ParasiticNode *node)
{
  return node_values_[node];
}

////////////////////////////////////////////////////////////////

class ReduceToPiElmore : public ReduceToPi
{
public:
  ReduceToPiElmore(StaState *sta);
  Parasitic *makePiElmore(const Parasitic *parasitic_network,
                          const Pin *drvr_pin,
                          ParasiticNode *drvr_node,
                          float coupling_cap_factor,
                          const RiseFall *rf,
                          const Scene *scene,
                          const MinMax *min_max);
  void reduceElmoreDfs(const Pin *drvr_pin,
                       ParasiticNode *node,
                       ParasiticResistor *from_res,
                       double elmore,
                       Parasitic *pi_elmore);
};

Parasitic *
reduceToPiElmore(const Parasitic *parasitic_network,
                 const Pin *drvr_pin,
                 const RiseFall *rf,
                 float coupling_cap_factor,
                 const Scene *scene,
                 const MinMax *min_max,
                 StaState *sta)
{
  Parasitics *parasitics = scene->parasitics(min_max);
  ParasiticNode *drvr_node =
    parasitics->findParasiticNode(parasitic_network, drvr_pin);
  if (drvr_node) {
    debugPrint(sta->debug(), "parasitic_reduce", 1, "Reduce driver {} {} {}",
               sta->network()->pathName(drvr_pin),
               rf->shortName(),
               min_max->to_string());
    ReduceToPiElmore reducer(sta);
    return reducer.makePiElmore(parasitic_network, drvr_pin, drvr_node,
                                coupling_cap_factor, rf, scene, min_max);
  }
  return nullptr;
}

ReduceToPiElmore::ReduceToPiElmore(StaState *sta) :
  ReduceToPi(sta)
{
}

Parasitic *
ReduceToPiElmore::makePiElmore(const Parasitic *parasitic_network,
                               const Pin *drvr_pin,
                               ParasiticNode *drvr_node,
                               float coupling_cap_factor,
                               const RiseFall *rf,
                               const Scene *scene,
                               const MinMax *min_max)
{
  float c2, rpi, c1;
  reduceToPi(parasitic_network, drvr_pin, drvr_node, coupling_cap_factor,
             rf, scene, min_max, c2, rpi, c1);
  Parasitic *pi_elmore = parasitics_->makePiElmore(drvr_pin, rf, min_max,
                                                   c2, rpi, c1);
  parasitics_->setIsReducedParasiticNetwork(pi_elmore, true);
  reduceElmoreDfs(drvr_pin, drvr_node, nullptr, 0.0, pi_elmore);
  return pi_elmore;
}

// Find elmore delays on 2nd DFS search using downstream capacitances
// set by reducePiDfs.
void
ReduceToPiElmore::reduceElmoreDfs(const Pin *drvr_pin,
                                  ParasiticNode *node,
                                  ParasiticResistor *from_res,
                                  double elmore,
                                  Parasitic *pi_elmore)
{
  const Pin *pin = parasitics_->pin(node);
  if (from_res && pin) {
    if (network_->isLoad(pin)) {
      debugPrint(debug_, "parasitic_reduce", 2, " Load {} elmore={:.3g}",
                 network_->pathName(pin),
                 elmore);
      parasitics_->setElmore(pi_elmore, pin, elmore);
    }
  }
  visit(node);
  ParasiticResistorSeq &resistors = resistor_map_[node];
  for (ParasiticResistor *resistor : resistors) {
    ParasiticNode *onode = parasitics_->otherNode(resistor, node);
    if (resistor != from_res
        && !isVisited(onode)
        && !isLoopResistor(resistor)) {
      float r = parasitics_->value(resistor);
      double onode_elmore = elmore + r * downstreamCap(onode);
      reduceElmoreDfs(drvr_pin, onode, resistor, onode_elmore, pi_elmore);
    }
  }
  leave(node);
}

////////////////////////////////////////////////////////////////

class ReduceToPiPoleResidue2 : public ReduceToPi
{
public:
  ReduceToPiPoleResidue2(StaState *sta);
  ~ReduceToPiPoleResidue2() override;
  void findPolesResidues(const Parasitic *parasitic_network,
                         Parasitic *pi_pole_residue,
                         const Pin *drvr_pin,
                         ParasiticNode *drvr_node);
  Parasitic *makePiPoleResidue2(const Parasitic *parasitic_network,
                                const Pin *drvr_pin,
                                ParasiticNode *drvr_node,
                                float coupling_cap_factor,
                                const RiseFall *rf,
                                const Scene *scene,
                                const MinMax *min_max);

private:
  void findMoments(const Pin *drvr_pin,
                   ParasiticNode *drvr_node,
                   int moment_count);
  void findMoments(const Pin *drvr_pin,
                   ParasiticNode *node,
                   double from_volt,
                   ParasiticResistor *from_res,
                   int moment_index);
  double findBranchCurrents(const Pin *drvr_pin,
                            ParasiticNode *node,
                            ParasiticResistor *from_res,
                            int moment_index);
  double moment(ParasiticNode *node,
                int moment_index);
  void setMoment(ParasiticNode *node,
                 double moment,
                 int moment_index);
  double current(ParasiticResistor *res);
  void setCurrent(ParasiticResistor *res,
                  double i);
  void findPolesResidues(Parasitic *pi_pole_residue,
                         const Pin *drvr_pin,
                         const Pin *load_pin,
                         ParasiticNode *load_node);

  // Resistor/capacitor currents.
  ResistorCurrentMap currents_;
  ParasiticNodeValueMap *moments_ {nullptr};
};

ReduceToPiPoleResidue2::ReduceToPiPoleResidue2(StaState *sta) :
  ReduceToPi(sta)
{
}

// The interconnect moments are found using RICE.
// "RICE: Rapid Interconnect Circuit Evaluation Using AWE",
// Curtis Ratzlaff and Lawrence Pillage, IEEE Transactions on
// Computer-Aided Design of Integrated Circuits and Systems,
// Vol 13, No 6, June 1994, pg 763-776.
//
// The poles and residues are found using these algorithms.
// "An Explicit RC-Circuit Delay Approximation Based on the First
// Three Moments of the Impulse Response", Proceedings of the 33rd
// Design Automation Conference, 1996, pg 611-616.
Parasitic *
reduceToPiPoleResidue2(const Parasitic *parasitic_network,
                       const Pin *drvr_pin,
                       const RiseFall *rf,
                       float coupling_cap_factor,
                       const Scene *scene,
                       const MinMax *min_max,
                       StaState *sta)
{
  Parasitics *parasitics = scene->parasitics(min_max);
  ParasiticNode *drvr_node =
    parasitics->findParasiticNode(parasitic_network, drvr_pin);
  if (drvr_node) {
    debugPrint(sta->debug(), "parasitic_reduce", 1, "Reduce driver {}",
               sta->network()->pathName(drvr_pin));
    ReduceToPiPoleResidue2 reducer(sta);
    return reducer.makePiPoleResidue2(parasitic_network, drvr_pin, drvr_node,
                                      coupling_cap_factor, rf,
                                      scene, min_max);
  }
  return nullptr;
}

Parasitic *
ReduceToPiPoleResidue2::makePiPoleResidue2(const Parasitic *parasitic_network,
                                           const Pin *drvr_pin,
                                           ParasiticNode *drvr_node,
                                           float coupling_cap_factor,
                                           const RiseFall *rf,
                                           const Scene *scene,
                                           const MinMax *min_max)
{
  float c2, rpi, c1;
  reduceToPi(parasitic_network, drvr_pin, drvr_node,
             coupling_cap_factor, rf, scene, min_max,
             c2, rpi, c1);
  Parasitic *pi_pole_residue = parasitics_->makePiPoleResidue(drvr_pin,
                                                              rf, min_max,
                                                              c2, rpi, c1);
  parasitics_->setIsReducedParasiticNetwork(pi_pole_residue, true);
  findPolesResidues(parasitic_network, pi_pole_residue, drvr_pin, drvr_node);
  return pi_pole_residue;
}

ReduceToPiPoleResidue2::~ReduceToPiPoleResidue2()
{
  delete [] moments_;
}

void
ReduceToPiPoleResidue2::findPolesResidues(const Parasitic *parasitic_network,
                                          Parasitic *pi_pole_residue,
                                          const Pin *drvr_pin,
                                          ParasiticNode *drvr_node)
{
  moments_ = new ParasiticNodeValueMap[4];
  findMoments(drvr_pin, drvr_node, 4);

  PinConnectedPinIterator *pin_iter = network_->connectedPinIterator(drvr_pin);
  while (pin_iter->hasNext()) {
    const Pin *pin = pin_iter->next();
    if (network_->isLoad(pin)) {
      ParasiticNode *load_node =
        parasitics_->findParasiticNode(parasitic_network, pin);
      if (load_node) {
        findPolesResidues(pi_pole_residue, drvr_pin, pin, load_node);
      }
    }
  }
  delete pin_iter;
}

void
ReduceToPiPoleResidue2::findMoments(const Pin *drvr_pin,
                                    ParasiticNode *drvr_node,
                                    int moment_count)
{
  // Driver model thevenin resistance.
  double rd = 0.0;
  // Zero'th moments are all 1 because Vin(0)=1 and there is no
  // current thru the resistors.  Thus, there is no point in doing a
  // pass to find the zero'th moments.
  for (int moment_index = 1; moment_index < moment_count; moment_index++) {
    double rd_i = findBranchCurrents(drvr_pin, drvr_node, nullptr, moment_index);
    double rd_volt = rd_i * rd;
    setMoment(drvr_node, 0.0, moment_index);
    findMoments(drvr_pin, drvr_node, -rd_volt, nullptr, moment_index);
  }
}

double
ReduceToPiPoleResidue2::findBranchCurrents(const Pin *drvr_pin,
                                           ParasiticNode *node,
                                           ParasiticResistor *from_res,
                                           int moment_index)
{
  visit(node);
  double branch_i = 0.0;
  double coupling_cap = 0.0;
  ParasiticResistorSeq &resistors = resistor_map_[node];
  for (ParasiticResistor *resistor : resistors) {
    ParasiticNode *onode = parasitics_->otherNode(resistor, node);
    // One commercial extractor creates resistors with identical from/to nodes.
    if (onode != node
        && resistor != from_res
        && !isVisited(onode)
        && !isLoopResistor(resistor)) {
      branch_i += findBranchCurrents(drvr_pin, onode, resistor, moment_index);
    }
  }
  ParasiticCapacitorSeq &capacitors = capacitor_map_[node];
  for (ParasiticCapacitor *capacitor : capacitors)
    coupling_cap += parasitics_->value(capacitor);

  double cap = parasitics_->nodeGndCap(node)
    + coupling_cap * coupling_cap_multiplier_
    + pinCapacitance(node);
  branch_i += cap * moment(node, moment_index - 1);
  leave(node);
  if (from_res) {
    setCurrent(from_res, branch_i);
    debugPrint(debug_, "parasitic_reduce", 3, " res i={:.3g}", branch_i);
  }
  return branch_i;
}

void
ReduceToPiPoleResidue2::findMoments(const Pin *drvr_pin,
                                    ParasiticNode *node,
                                    double from_volt,
                                    ParasiticResistor *from_res,
                                    int moment_index)
{
  visit(node);
  ParasiticResistorSeq &resistors = resistor_map_[node];
  for (ParasiticResistor *resistor : resistors) {
    ParasiticNode *onode = parasitics_->otherNode(resistor, node);
    // One commercial extractor creates resistors with identical from/to nodes.
    if (onode != node
        && resistor != from_res
        && !isVisited(onode)
        && !isLoopResistor(resistor)) {
      double r = parasitics_->value(resistor);
      double r_volt = r * current(resistor);
      double onode_volt = from_volt - r_volt;
      setMoment(onode, onode_volt, moment_index);
      debugPrint(debug_, "parasitic_reduce", 3, " moment {} {} {:.3g}",
                 parasitics_->name(onode),
                 moment_index,
                 onode_volt);
      findMoments(drvr_pin, onode, onode_volt, resistor, moment_index);
    }
  }
  leave(node);
}

double
ReduceToPiPoleResidue2::moment(ParasiticNode *node,
                               int moment_index)
{
  // Zero'th moments are all 1.
  if (moment_index == 0)
    return 1.0;
  else {
    ParasiticNodeValueMap &map = moments_[moment_index];
    return map[node];
  }
}

void
ReduceToPiPoleResidue2::setMoment(ParasiticNode *node,
                                  double moment,
                                  int moment_index)
{
  // Zero'th moments are all 1.
  if (moment_index > 0) {
    ParasiticNodeValueMap &map = moments_[moment_index];
    map[node] = moment;
  }
}

double
ReduceToPiPoleResidue2::current(ParasiticResistor *res)
{
  return currents_[res];
}

void
ReduceToPiPoleResidue2::setCurrent(ParasiticResistor *res,
                                   double i)
{
  currents_[res] = i;
}

void
ReduceToPiPoleResidue2::findPolesResidues(Parasitic *pi_pole_residue,
                                          const Pin *,
                                          const Pin *load_pin,
                                          ParasiticNode *load_node)
{
  double m1 = moment(load_node, 1);
  double m2 = moment(load_node, 2);
  double m3 = moment(load_node, 3);
  double p1 = -m2 / m3;
  double p2 = p1 * (1.0 /  m1 - m1 / m2) / (m1 / m2 - m2 / m3);
  if (p1 <= 0.0
      || p2 <= 0.0
      // Coincident poles.  Not handled by delay calculator.
      || p1 == p2
      || m1 / m2 == m2 / m3) {
    double p1 = -1.0 / m1;
    double k1 = 1.0;
    debugPrint(debug_, "parasitic_reduce", 3, " load {} p1={:.3g} k1={:.3g}",
               network_->pathName(load_pin), p1, k1);
    ComplexFloatSeq *poles = new ComplexFloatSeq(1);
    ComplexFloatSeq *residues = new ComplexFloatSeq(1);
    (*poles)[0] = ComplexFloat(p1, 0.0);
    (*residues)[0] = ComplexFloat(k1, 0.0);
    parasitics_->setPoleResidue(pi_pole_residue, load_pin, poles, residues);
  }
  else {
    double k1 = p1 * p1 * (1.0 + m1 * p2) / (p1 - p2);
    double k2 = -p2 * p2 * (1.0 + m1 * p1) / (p1 - p2);
    if (k1 < 0.0 && k2 > 0.0) {
      // Swap p1 and p2.
      double p = p2, k = k2;
      p2 = p1;
      k2 = k1;
      p1 = p;
      k1 = k;
    }
    debugPrint(debug_, "parasitic_reduce", 3,
               " load {} p1={:.3g} p2={:.3g} k1={:.3g} k2={:.3g}",
               network_->pathName(load_pin), p1, p2, k1, k2);

    ComplexFloatSeq *poles = new ComplexFloatSeq(2);
    ComplexFloatSeq *residues = new ComplexFloatSeq(2);
    (*poles)[0] = ComplexFloat(p1, 0.0);
    (*residues)[0] = ComplexFloat(k1, 0.0);
    (*poles)[1] = ComplexFloat(p2, 0.0);
    (*residues)[1] = ComplexFloat(k2, 0.0);
    parasitics_->setPoleResidue(pi_pole_residue, load_pin, poles, residues);
  }
}

} // namespace sta
