text:
fluidData["N2"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "N2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7727-37-9",
                       meltingPoint =            63.15,
                       normalBoilingPoint =      77.35,
                       criticalTemperature =    126.20,
                       criticalPressure =        33.98e5,
                       criticalMolarVolume =     90.10e-6,
                       acentricFactor =           0.037,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["O2"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "O2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7782-44-7",
                       meltingPoint =            54.36,
                       normalBoilingPoint =      90.17,
                       criticalTemperature =    154.58,
                       criticalPressure =        50.43e5,
                       criticalMolarVolume =     73.37e-6,
                       acentricFactor =         0.022,
                       dipoleMoment =           0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["CL2"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "CL2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7782-50-5",
                       meltingPoint =           172.19,
                       normalBoilingPoint =     239.12,
                       criticalTemperature =    417.00,
                       criticalPressure =        77.00e5,
                       criticalMolarVolume =    124.00e-6,
                       acentricFactor =         0.069,
                       dipoleMoment =           0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["F2"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "F2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7782-41-4",
                       meltingPoint =            53.48,
                       normalBoilingPoint =      84.95,
                       criticalTemperature =    144.30,
                       criticalPressure =        52.15e5,
                       criticalMolarVolume =     66.20e-6,
                       acentricFactor =           0.051,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["CO2"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "CO2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "124-38-9",
                       meltingPoint =           216.58,
                       normalBoilingPoint =     -1.0,
                       criticalTemperature =    304.12,
                       criticalPressure =        73.74e5,
                       criticalMolarVolume =     94.07e-6,
                       acentricFactor =           0.225,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["CO"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "CO",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "630-08-0",
                       meltingPoint =            68.15,
                       normalBoilingPoint =      81.66,
                       criticalTemperature =    132.85,
                       criticalPressure =        34.94e5,
                       criticalMolarVolume =     93.10e-6,
                       acentricFactor =           0.045,
                       dipoleMoment =             0.1,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["H2"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "H2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "800000-51-5",
                       meltingPoint =            13.56,
                       normalBoilingPoint =      20.38,
                       criticalTemperature =     33.25,
                       criticalPressure =        12.97e5,
                       criticalMolarVolume =     65.00e-6,
                       acentricFactor =          -0.216,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["H2O"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "H2O",
                       iupacName =              "oxidane",
                       structureFormula =       "H2O",
                       casRegistryNumber =      "7732-18-5",
                       meltingPoint =           273.15,
                       normalBoilingPoint =     373.124,
                       criticalTemperature =    647.096,
                       criticalPressure =       220.64e5,
                       criticalMolarVolume =     55.95e-6,
                       acentricFactor =           0.344,
                       dipoleMoment =             1.8,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["N2O"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "N2O",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "10024-97-2",
                       meltingPoint =           182.33,
                       normalBoilingPoint =     184.67,
                       criticalTemperature =    309.60,
                       criticalPressure =        72.55e5,
                       criticalMolarVolume =     97.00e-6,
                       acentricFactor =         0.142,
                       dipoleMoment =           0.2,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["NO"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "NO",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "10102-43-9",
                       meltingPoint =           109.51,
                       normalBoilingPoint =     121.38,
                       criticalTemperature =    180.00,
                       criticalPressure =        64.80e5,
                       criticalMolarVolume =     58.00e-6,
                       acentricFactor =           0.582,
                       dipoleMoment =             0.2,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["NO2"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "NO2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "10102-44-0",
                       meltingPoint =           261.95,
                       normalBoilingPoint =     294.0,
                       criticalTemperature =    431.35,
                       criticalPressure =       101.33e5,
                       criticalMolarVolume =    82.5e-6,
                       acentricFactor =         0.849,
                       dipoleMoment =           0.32,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["NH3"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "NH3",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7664-41-7",
                       meltingPoint =           195.41,
                       normalBoilingPoint =     239.82,
                       criticalTemperature =    405.40,
                       criticalPressure =       113.53e5,
                       criticalMolarVolume =     72.47e-6,
                       acentricFactor =           0.257,
                       dipoleMoment =             1.5,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["SO2"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "SO2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7446-09-5",
                       meltingPoint =           197.67,
                       normalBoilingPoint =     263.13,
                       criticalTemperature =    430.80,
                       criticalPressure =        78.84e5,
                       criticalMolarVolume =    122.00e-6,
                       acentricFactor =         0.245,
                       dipoleMoment =           1.6,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["SO3"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "SO3",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7446-11-9",
                       meltingPoint =           289.95,
                       normalBoilingPoint =     317.90,
                       criticalTemperature =    490.90,
                       criticalPressure =        82.10e5,
                       criticalMolarVolume =    126.50e-6,
                       acentricFactor =         0.422,
                       dipoleMoment =           0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["Ar"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "Ar",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7440-37-1",
                       meltingPoint =            83.80,
                       normalBoilingPoint =      87.27,
                       criticalTemperature =    150.86,
                       criticalPressure =        48.98e5,
                       criticalMolarVolume =     74.57e-6,
                       acentricFactor =          -0.002,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["He"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "He",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7440-59-7",
                       meltingPoint =             2.15,
                       normalBoilingPoint =       4.30,
                       criticalTemperature =      5.19,
                       criticalPressure =         2.27e5,
                       criticalMolarVolume =     57.30e-6,
                       acentricFactor =          -0.390,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["Ne"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "Ne",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "7440-01-9",
                       meltingPoint =            24.56,
                       normalBoilingPoint =      27.07,
                       criticalTemperature =     44.40,
                       criticalPressure =        27.60e5,
                       criticalMolarVolume =     41.70e-6,
                       acentricFactor =          -0.016,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["CH4"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "CH4",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "74-82-8",
                       meltingPoint =            90.69,
                       normalBoilingPoint =     111.66,
                       criticalTemperature =    190.56,
                       criticalPressure =        45.99e5,
                       criticalMolarVolume =     98.60e-6,
                       acentricFactor =           0.011,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C2H6"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C2H6",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "74-84-0",
                       meltingPoint =            90.35,
                       normalBoilingPoint =     184.55,
                       criticalTemperature =    305.32,
                       criticalPressure =        48.72e5,
                       criticalMolarVolume =    145.50e-6,
                       acentricFactor =           0.099,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C3H8"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C3H8",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "74-98-6",
                       meltingPoint =            91.45,
                       normalBoilingPoint =     231.02,
                       criticalTemperature =    369.83,
                       criticalPressure =        42.48e5,
                       criticalMolarVolume =    200.00e-6,
                       acentricFactor =           0.152,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C4H10_n_butane"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C4H10",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "106-97-8",
                       meltingPoint =           134.79,
                       normalBoilingPoint =     272.66,
                       criticalTemperature =    425.12,
                       criticalPressure =        37.96e5,
                       criticalMolarVolume =    255.00e-6,
                       acentricFactor =           0.20,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C5H12_n_pentane"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C5H12",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "109-66-0",
                       meltingPoint =           143.43,
                       normalBoilingPoint =     309.22,
                       criticalTemperature =    469.70,
                       criticalPressure =        33.70e5,
                       criticalMolarVolume =    311.00e-6,
                       acentricFactor =           0.252,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C6H14_n_hexane"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C6H14",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "110-54-3",
                       meltingPoint =           177.84,
                       normalBoilingPoint =     341.88,
                       criticalTemperature =    507.60,
                       criticalPressure =        30.25e5,
                       criticalMolarVolume =    368.00e-6,
                       acentricFactor =           0.300,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C7H16_n_heptane"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C7H16",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "142-82-5",
                       meltingPoint =           182.59,
                       normalBoilingPoint =     371.57,
                       criticalTemperature =    540.20,
                       criticalPressure =        27.40e5,
                       criticalMolarVolume =    428.00e-6,
                       acentricFactor =           0.350,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C2H4"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C2H4",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "74-85-1",
                       meltingPoint =           103.99,
                       normalBoilingPoint =     169.42,
                       criticalTemperature =    282.34,
                       criticalPressure =        50.41e5,
                       criticalMolarVolume =    131.10e-6,
                       acentricFactor =           0.087,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C3H6_propylene"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C3H6",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "115-07-1",
                       meltingPoint =            87.89,
                       normalBoilingPoint =     225.46,
                       criticalTemperature =    364.90,
                       criticalPressure =        46.00e5,
                       criticalMolarVolume =    184.60e-6,
                       acentricFactor =           0.142,
                       dipoleMoment =             0.4,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C4H8_1_butene"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C4H8",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "106-98-9",
                       meltingPoint =            87.79,
                       normalBoilingPoint =     266.92,
                       criticalTemperature =    419.50,
                       criticalPressure =        40.20e5,
                       criticalMolarVolume =    240.80e-6,
                       acentricFactor =           0.194,
                       dipoleMoment =             0.3,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C5H10_1_pentene"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C5H10",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "109-67-1",
                       meltingPoint =           106.95,
                       normalBoilingPoint =     303.11,
                       criticalTemperature =    464.80,
                       criticalPressure =        35.60e5,
                       criticalMolarVolume =    298.40e-6,
                       acentricFactor =           0.237,
                       dipoleMoment =             0.4,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C6H12_1_hexene"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C6H12",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "592-41-6",
                       meltingPoint =           133.34,
                       normalBoilingPoint =     336.63,
                       criticalTemperature =    504.00,
                       criticalPressure =        31.43e5,
                       criticalMolarVolume =    355.10e-6,
                       acentricFactor =           0.281,
                       dipoleMoment =             0.4,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C7H14_1_heptene"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C7H14",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "592-76-7",
                       meltingPoint =           153.45,
                       normalBoilingPoint =     366.79,
                       criticalTemperature =    537.30,
                       criticalPressure =        29.20e5,
                       criticalMolarVolume =    409.00e-6,
                       acentricFactor =           0.343,
                       dipoleMoment =             0.3,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C2H2_vinylidene"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C2H2",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "74-86-2",
                       meltingPoint =           192.35,
                       normalBoilingPoint =     188.40,
                       criticalTemperature =    308.30,
                       criticalPressure =        61.14e5,
                       criticalMolarVolume =    112.20e-6,
                       acentricFactor =           0.189,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C6H6"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C6H6",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "71-43-2",
                       meltingPoint =           278.68,
                       normalBoilingPoint =     353.24,
                       criticalTemperature =    562.05,
                       criticalPressure =        48.95e5,
                       criticalMolarVolume =    256.00e-6,
                       acentricFactor =           0.210,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C8H18_n_octane"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C8H18",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "111-65-9",
                       meltingPoint =           216.39,
                       normalBoilingPoint =     398.82,
                       criticalTemperature =    568.70,
                       criticalPressure =        24.90e5,
                       criticalMolarVolume =    492.00e-6,
                       acentricFactor =           0.399,
                       dipoleMoment =             0.0,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C8H10_ethylbenz"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C8H10",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "100-41-4",
                       meltingPoint =           178.18,
                       normalBoilingPoint =     409.36,
                       criticalTemperature =    617.15,
                       criticalPressure =        36.09e5,
                       criticalMolarVolume =    374.00e-6,
                       acentricFactor =           0.304,
                       dipoleMoment =             0.4,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["CH3OH"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "CH3OH",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "67-56-1",
                       meltingPoint =           175.49,
                       normalBoilingPoint =     337.69,
                       criticalTemperature =    512.64,
                       criticalPressure =        80.97e5,
                       criticalMolarVolume =    118.00e-6,
                       acentricFactor =           0.565,
                       dipoleMoment =             1.7,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C2H5OH"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C2H5OH",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "64-17-5",
                       meltingPoint =           159.05,
                       normalBoilingPoint =     351.80,
                       criticalTemperature =    513.92,
                       criticalPressure =        61.48e5,
                       criticalMolarVolume =    167.00e-6,
                       acentricFactor =           0.649,
                       dipoleMoment =             1.7,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

#= no corresponding data in Modelica.Media.IdealGases.Common.SingleGasesData.jl
fluidData["C3H7OH"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C3H7OH",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "71-23-8",
                       meltingPoint =           147.00,
                       normalBoilingPoint =     370.93,
                       criticalTemperature =    536.78,
                       criticalPressure =        51.75e5,
                       criticalMolarVolume =    219.00e-6,
                       acentricFactor =           0.629,
                       dipoleMoment =             1.7,
                       molarMass =              60.1e-3,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);

fluidData["C4H9OH"] = ModiaMedia.IdealGasFluidConstants(
                       chemicalFormula =        "C4H9OH",
                       iupacName =              "unknown",
                       structureFormula =       "unknown",
                       casRegistryNumber =      "71-36-3",
                       meltingPoint =           183.35,
                       normalBoilingPoint =     390.88,
                       criticalTemperature =    563.05,
                       criticalPressure =        44.23e5,
                       criticalMolarVolume =    275.00e-6,
                       acentricFactor =           0.589,
                       dipoleMoment =             1.8,
                       molarMass =              74.12e-3,
                       hasDipoleMoment =       true,
                       hasIdealGasHeatCapacity=true,
                       hasCriticalData =       true,
                       hasAcentricFactor =     true);
=#


text:
@testset "Floyd Warshall" begin
    g3 = PathGraph(5)
    d = [0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]
    for g in testgraphs(g3)
      z = @inferred(floyd_warshall_shortest_paths(g, d))
      @test z.dists[3, :][:] == [7, 6, 0, 11, 27]
      @test z.parents[3, :][:] == [2, 3, 0, 3, 4]

      @test @inferred(enumerate_paths(z))[2][2] == []
      @test @inferred(enumerate_paths(z))[2][4] == enumerate_paths(z, 2)[4] == enumerate_paths(z, 2, 4) == [2, 3, 4]
    end
    g4 = PathDiGraph(4)
    d = ones(4, 4)
    for g in testdigraphs(g4)
      z = @inferred(floyd_warshall_shortest_paths(g, d))
      @test length(enumerate_paths(z, 4, 3)) == 0
      @test length(enumerate_paths(z, 4, 1)) == 0
      @test length(enumerate_paths(z, 2, 3)) == 2
    end
end

ext: jl


text:
## mlutils.jl - a set of ML utilities for parameter initalization, performance metrics, etc.

## initialization

function xavier_init(dim_in, dim_out; c=1)
    low = -c * sqrt(6.0 / (dim_in + dim_out))
    high = c * sqrt(6.0 / (dim_in + dim_out))
    return rand(Uniform(low, high), dim_in, dim_out)
end

## one-hot encoding

function class_to_index(classes)
    classes = sort(classes)
    return Dict(c => i for (i, c) in enumerate(classes))
end

function one_hot(::Type{T}, num_classes::Int, vals::AbstractVector{<:Integer}) where T
    ret = zeros(T, num_classes, length(vals))
    for (i, v) in enumerate(vals)
        ret[v, i] = T(1)
    end
    return ret
end

one_hot(num_classes::Integer, vals) = one_hot(Float64, num_classes, vals)

function one_hot(::Type{T}, c2i::Dict, vals) where T
    result = zeros(T, length(c2i), length(vals))
    for (j, val) in enumerate(vals)
        i = c2i[val]
        result[i, j] = one(T)
    end
    return result
end

one_hot(c2i::Dict, vals) = one_hot(Float64, c2i, vals)

## metrics

function accuracy(ŷ::AbstractMatrix, y::AbstractMatrix)
    return mean(findmax(@view ŷ[:,i])[2] == findmax(@view y[:,i])[2] for i=1:size(ŷ,2))
end

## CUDA support

@require CuArrays begin

    """
    Convert to a corresponding CUDA object.

    For arrays returns an instance of CuArray.
    For structs returns similar struct with all arrays fields converted to CuArray.
    """
    to_cuda(a::AbstractArray{T,N}) where {T,N} = convert(CuArray{Float32,N}, a)
    to_cuda(a::CuArray) = a

    function to_cuda(m)
        if isstruct(m)
            new_m = Espresso.struct_like(m)
            for fld in fieldnames(typeof(m))
                v = getfield(m, fld)
                new_v = to_cuda(v)
                setfield!(new_m, fld, new_v)
            end
            return new_m
        else
            return m
        end
    end

end


text:
"""
$(TYPEDEF)

Defines a time dependent Hamiltonian object with sparse Matrices.

# Fields

$(FIELDS)
"""
struct SparseHamiltonian{T<:Number} <: AbstractSparseHamiltonian{T}
    "List of time dependent functions"
    f::Any
    "List of constant matrices"
    m::Vector{SparseMatrixCSC{T,Int}}
    "Internal cache"
    u_cache::SparseMatrixCSC{T,Int}
    "Size"
    size::Tuple
    "Eigen decomposition routine"
    EIGS::Any
end

"""
$(SIGNATURES)

Constructor of the `SparseHamiltonian` type. `funcs` and `mats` are lists of time-dependent functions and the corresponding matrices. The Hamiltonian can be represented as ``∑ᵢfuncs[i](s)×mats[i]``. `unit` specifies wether `:h` or `:ħ` is set to one when defining `funcs` and `mats`. The `mats` will be scaled by ``2π`` if unit is `:h`.

`EIGS` is the initializer for the eigendecomposition routine for the Hamiltonian. It should return a function of the signature: `(H, s, lvl) -> (w, v)`. The default method `EIGEN_DEFAULT` will use `LAPACK` routine.
"""
function SparseHamiltonian(funcs, mats; unit = :h, EIGS = EIGEN_DEFAULT)
    if any((x) -> size(x) != size(mats[1]), mats)
        throw(ArgumentError("Matrices in the list do not have the same size."))
    end
    if is_complex(funcs, mats)
        mats = complex.(mats)
    end
    cache = similar(sum(mats))
    fill!(cache, 0.0)
    mats = unit_scale(unit) * mats
    EIGS = EIGS(cache)
    SparseHamiltonian(funcs, mats, cache, size(mats[1]), EIGS)
end

"""
    function (h::SparseHamiltonian)(t::Real)

Calling the Hamiltonian returns the value ``2πH(t)``.
"""
function (h::SparseHamiltonian)(t::Real)
    fill!(h.u_cache, 0.0)
    for (f, m) in zip(h.f, h.m)
        h.u_cache .+= f(t) * m
    end
    h.u_cache
end

# The argument `p` is not essential for `SparseHamiltonian`
# It exists to keep the `update_cache!` interface consistent across
# all `AbstractHamiltonian` types
function update_cache!(cache, H::SparseHamiltonian, tf, s::Real)
    fill!(cache, 0.0)
    for (f, m) in zip(H.f, H.m)
        cache .+= -1.0im * f(s) * m
    end
end

function update_vectorized_cache!(cache, H::SparseHamiltonian, tf, s::Real)
    hmat = H(s)
    iden = sparse(I, size(H))
    cache .= 1.0im * (transpose(hmat) ⊗ iden - iden ⊗ hmat)
end

function (h::SparseHamiltonian)(
    du,
    u,
    p,
    s::Real,
) where {T<:Number}
    H = h(s)
    du .= -1.0im * (H * u - u * H)
end

function Base.convert(S::Type{T}, H::SparseHamiltonian{M}) where {T<:Complex, M}
    mats = [convert.(S, x) for x in H.m]
    cache = similar(H.u_cache, complex{M})
    SparseHamiltonian(H.f, mats, cache, size(H))
end

function Base.convert(S::Type{T}, H::SparseHamiltonian{M}) where {T<:Real, M}
    f_val = sum((x) -> x(0.0), H.f)
    if !(typeof(f_val) <: Real)
        throw(TypeError(:convert, "H.f", Real, typeof(f_val)))
    end
    mats = [convert.(S, x) for x in H.m]
    cache = similar(H.u_cache, real(M))
    SparseHamiltonian(H.f, mats, cache, size(H), H.EIGS)
end


text:
export assembleStiffAndForce,assembleMassMatrix!,assembleInternalForce, getEqns

@doc raw"""
    assembleInternalForce(domain::Domain, Δt::Float64 = 0.0)

Computes the internal force vector $F_\mathrm{int}$ of length `neqs`
- `globdat`: GlobalData
- `domain`: Domain, finite element domain, for data structure
- `Δt`:  Float64, current time step size

Only the information in `domain` is used for computing internal force.
Therefore, the boundary conditions in `domain` must be set appropriately.
"""
function assembleInternalForce(domain::Domain, Δt::Float64 = 0.0)
    Fint = zeros(Float64, domain.neqs)
    neles = domain.neles

    # Loop over the elements in the elementGroup
    for iele  = 1:neles
        element = domain.elements[iele]

      # Get the element nodes
        el_nodes = getNodes(element)

      # Get the element nodes
        el_eqns = getEqns(domain, iele)

        el_dofs = getDofs(domain, iele)

        el_state  = getState(domain, el_dofs)

        el_Dstate = getDstate(domain, el_dofs)

      # Get the element contribution by calling the specified action
        fint = getInternalForce(element, el_state, el_Dstate, Δt)

      # Assemble in the global array
        el_eqns_active = (el_eqns .>= 1)
        Fint[el_eqns[el_eqns_active]] += fint[el_eqns_active]
    end

    return Fint
end
assembleInternalForce(globdat::GlobalData, domain::Domain, Δt::Float64 = 0.0) = assembleInternalForce(domain, Δt)

@doc raw"""
    assembleInternalForce(domain::Domain, nn::Function, E_all::PyObject, DE_all::PyObject, w∂E∂u_all::PyObject, σ0_all::PyObject)

Computes local internal force f_int and then assemble to F_int, which generates inverse problem automatically.

- `domain`: finite element domain
- `nn`: constitutive relation for expressing `stress = f(strain)`, assuming `stress` and `strain` are defined on Gauss points (`(neles*nGauss) × nstrains`).
- `E_all`: strain data of size `(neles*nGauss) × nstrains` at the **current** time step.
- `DE_all`: strain data of size `(neles*nGauss) × nstrains` at the **last** time step.
- `w∂E∂u_all`: sensitivity matrix of size `(neles*nGauss) x ndofs_per_element x nstrains`; `neles*nGauss` is the number of Gaussian quadrature points,
  `ndofs_per_element` is the number of freedoms per element, and `nstrain` is the number of strain components.
  The sensitivity matrix already considers the quadrature weights.

```math
s_{g,j,i}^e = w_g^e\frac{\partial \epsilon_g^e}{\partial u_j^e}
```
where $e$ is the element index, $g$ is the Gaussian index.

- `σ0_all`: stress data of size `neles*nGauss×nstrains` at the **last** time step.

Return:

- $F_{\mathrm{int}}$:  internal force vector of length `neqns`
- $\sigma_{\mathrm{all}}$: predicted stress at **current** step, a matrix of size `(neles*nGauss) × nstrains`
"""
function assembleInternalForce(domain::Domain, nn::Function,
        E_all::Union{Array{Float64, 2}, PyObject}, DE_all::Union{Array{Float64, 2}, PyObject},
        w∂E∂u_all::Union{Array{Float64, 3}, PyObject}, σ0_all::Union{Array{Float64, 2}, PyObject})

    E_all, DE_all, w∂E∂u_all, σ0_all = convert_to_tensor([E_all, DE_all, w∂E∂u_all, σ0_all], [Float64, Float64, Float64, Float64])
    neles = domain.neles
    nGauss = length(domain.elements[1].weights)
    neqns_per_elem = length(getEqns(domain, 1))
    nstrains = size(E_all, 2)

    @assert size(E_all) == (neles * nGauss, nstrains)
    @assert size(DE_all) == (neles * nGauss, nstrains)
    @assert size(σ0_all) == (neles * nGauss, nstrains)

  # el_eqns_all, the equation numbers related to the Gaussian points in the element, negative value means Drichlet boundary
    el_eqns_all = zeros(Int32, neles * nGauss, neqns_per_elem)
  # el_eqns_active_all = el_eqns_all > 0
    el_eqns_active_all = zeros(Bool, neles * nGauss, neqns_per_elem)

  # Loop over the elements in the elementGroup to construct el_eqns_active_all , el_eqns_all and w∂E∂u_all
    for iele  = 1:neles
    # Get the element nodes
        el_eqns = getEqns(domain, iele)

    # Assemble in the global array
        el_eqns_active_all[(iele - 1) * nGauss + 1:iele * nGauss,:] = repeat((el_eqns .>= 1)', nGauss, 1)
        el_eqns_all[(iele - 1) * nGauss + 1:iele * nGauss,:] = repeat(el_eqns', nGauss, 1)
    end

  # get stress at each Gaussian points
    σ_all = nn(E_all, DE_all, σ0_all)

  # compute fint at each Gaussian quadrature points, fints[igp] = w∂E∂u_all[igp] * σ_all[igp]
    fints = squeeze(tf.matmul(w∂E∂u_all, tf.expand_dims(σ_all, 2)))

  # call the Cpp assembler to construct Fint
    cpp_fint = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib", "fint_comp")
    Fint = cpp_fint(fints, constant(el_eqns_all, dtype = Int32), constant(domain.neqs, dtype = Int32))

    return Fint, σ_all
end

@doc raw"""
    assembleStiffAndForce(domain::Domain, Δt::Float64 = 0.0)

Computes the internal force and stiffness matrix.

- `domain`: Domain, finite element domain, for data structure
- `Δt`:  Float64, current time step size

Returns a length `neqs` vector $F_{\mathrm{int}}$ and `neqs×neqs` sparse stiffness matrix.
"""
function assembleStiffAndForce(domain::Domain, Δt::Float64 = 0.0)
    Fint = zeros(Float64, domain.neqs)
  # K = zeros(Float64, domain.neqs, domain.neqs)
    ii = Int64[]; jj = Int64[]; vv = Float64[]
    neles = domain.neles

  # Loop over the elements in the elementGroup
    for iele  = 1:neles
        element = domain.elements[iele]

    # Get the element nodes
        el_nodes = getNodes(element)

    # Get the element nodes equation numbers
        el_eqns = getEqns(domain, iele)

    # Get the element nodes dof numbers
        el_dofs = getDofs(domain, iele)

        el_state  = getState(domain, el_dofs)

        el_Dstate = getDstate(domain, el_dofs)

    # Get the element contribution by calling the specified action
        fint, stiff  = getStiffAndForce(element, el_state, el_Dstate, Δt)

    # Assemble in the global array
        el_eqns_active = el_eqns .>= 1
    # K[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += stiff[el_eqns_active,el_eqns_active]
        Slocal = stiff[el_eqns_active,el_eqns_active]
        Idx = el_eqns[el_eqns_active]
        for i = 1:length(Idx)
            for j = 1:length(Idx)
                push!(ii, Idx[i])
                push!(jj, Idx[j])
                push!(vv, Slocal[i,j])
            end
        end
        Fint[el_eqns[el_eqns_active]] += fint[el_eqns_active]
    end
    Ksparse = sparse(ii, jj, vv, domain.neqs, domain.neqs)

    return Fint, Ksparse
end

assembleStiffAndForce(globdat::GlobalData, domain::Domain, Δt::Float64 = 0.0) = assembleStiffAndForce(domain, Δt)

@doc raw"""
    assembleMassMatrix!(globaldat::GlobalData, domain::Domain)

Computes the constant sparse mass matrix $M_{\mathrm{mass}}$, the lumped mass matrix $M_{\mathrm{lump}}$
due to time-dependent Dirichlet boundary conditions, and store them in `globaldat`.

```math
M_{\mathrm{mass}}\begin{bmatrix}
M & M_{ID}\\
M_{ID}^T & M_{DD}
\end{bmatrix}
```

Here M is a `neqns×neqns` matrix, and $M_{ID}$ is a `neqns×nd` matrix. $M_{\mathrm{lump}}$ assumes that the local mass matrix is a diagonal matrix.

- `globdat`: `GlobalData`
- `domain`: `Domain`, finite element domain, for data structure

![](./assets/massmatrix.png)
"""
function assembleMassMatrix!(globaldat::GlobalData, domain::Domain)
    if size(globaldat.M,1)>0
        return
    end
    Mlumped = zeros(Float64, domain.neqs)
    # M = zeros(Float64, domain.neqs, domain.neqs)
    iiM = Int64[]; jjM = Int64[]; vvM = Float64[]
    Mlumped = zeros(Float64, domain.neqs)
    neles = domain.neles

    nnodes, ndims = domain.nnodes, domain.ndims

    # MID = zeros(domain.neqs, sum(domain.EBC .== -2))
    iiMID = Int64[]; jjMID = Int64[]; vvMID = Float64[]

    # construct map from freedoms(first direction, second direction ...) to time-dependent Dirichlet freedoms
    dofs_to_EBCdofs = zeros(Int64, nnodes*ndims)
    if globaldat.EBC_func != nothing
        dof_id = 0
        for idof = 1:ndims
            for inode = 1:nnodes
                if (domain.EBC[inode, idof] == -2)
                    dof_id += 1
                    dofs_to_EBCdofs[inode + (idof - 1) * nnodes] = dof_id
                end
            end
        end
    end

    # Loop over the elements in the elementGroup
    for iele = 1:neles
        element = domain.elements[iele]

        # Get the element nodes
        el_nodes = getNodes(element)

        # Get the element nodes equation numbers
        el_eqns = getEqns(domain, iele)

        # Get the element nodes dof numbers
        el_dofs = getDofs(domain, iele)

        # Get the element contribution by calling the specified action
        lM, lumped = getMassMatrix(element)

        # Assemble in the global array

        el_eqns_active = (el_eqns .>= 1)

        # time-dependent Dirichlet boundary condition
        el_eqns_acc_active = (el_eqns .== -2)

        # M[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += lM[el_eqns_active, el_eqns_active]
        Idx = el_eqns[el_eqns_active]
        Mlocal = lM[el_eqns_active, el_eqns_active]
        for i = 1:length(Idx)
            for j = 1:length(Idx)
                push!(iiM, Idx[i])
                push!(jjM, Idx[j])
                push!(vvM, Mlocal[i,j])
            end
        end

        Mlumped[el_eqns[el_eqns_active]] += lumped[el_eqns_active]

        if globaldat.EBC_func != nothing
            Idx = el_eqns[el_eqns_active]
            Idy = dofs_to_EBCdofs[el_dofs[el_eqns_acc_active]]

            Mlocal = lM[el_eqns_active, el_eqns_acc_active]

            for i = 1:length(Idx)
                for j = 1:length(Idy)
                    push!(iiMID, Idx[i])
                    push!(jjMID, Idy[j])
                    push!(vvMID, Mlocal[i,j])
                end
            end
        end

    end

    globaldat.M = sparse(iiM, jjM, vvM, domain.neqs, domain.neqs)
    globaldat.Mlumped = Mlumped
    globaldat.MID = sparse(iiMID, jjMID, vvMID, domain.neqs, sum(domain.EBC .== -2))

end


text:
module ParametricOptInterface

using MathOptInterface
using DataStructures: OrderedDict

const MOI = MathOptInterface

const PARAMETER_INDEX_THRESHOLD = 1_000_000_000_000_000_000
const SUPPORTED_SETS =
    (MOI.LessThan{Float64}, MOI.EqualTo{Float64}, MOI.GreaterThan{Float64})
const SUPPORTED_VECTOR_SETS = (MOI.Nonnegatives, MOI.SecondOrderCone)

"""
    Parameter(::Float64)

The `Parameter` structure stores the numerical value associated to a given
parameter.

## Example

```julia-repl
julia> ParametricOptInterface.Parameter(5)
ParametricOptInterface.Parameter(5)
```
"""
struct Parameter <: MOI.AbstractScalarSet
    val::Float64
end

"""
    Optimizer{T, OT <: MOI.ModelLike} <: MOI.AbstractOptimizer

Declares a `Optimizer`, which allows the handling of parameters in a
optimization model.

## Example

```julia-repl
julia> ParametricOptInterface.Optimizer(GLPK.Optimizer())
ParametricOptInterface.Optimizer{Float64,GLPK.Optimizer}
```
"""
mutable struct Optimizer{T,OT<:MOI.ModelLike} <: MOI.AbstractOptimizer
    optimizer::OT
    parameters::Dict{MOI.VariableIndex,T}
    parameters_name::Dict{MOI.VariableIndex,String}
    updated_parameters::Dict{MOI.VariableIndex,T}
    variables::Dict{MOI.VariableIndex,MOI.VariableIndex}
    last_variable_index_added::Int64
    last_parameter_index_added::Int64
    # Store reference to affine constraints with parameters: v + p >= 0
    affine_constraint_cache::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.ScalarAffineTerm{Float64}},
    }
    # Store reference to parameter * variable constraints: p*v >= 0.0
    quadratic_constraint_cache_pv::Dict{
        MOI.ConstraintIndex,
        Vector{MOI.ScalarQuadraticTerm{Float64}},
    }
    # Store reference to parameter * parameter constraints: p*p >= var
    quadratic_constraint_cache_pp::Dict{
        MOI.ConstraintIndex,
        Vector{MOI.ScalarQuadraticTerm{Float64}},
    }
    # Store reference to constraints quad_variable_term + affine_with_parameters: v*v + p >= 0
    quadratic_constraint_cache_pc::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.ScalarAffineTerm{Float64}},
    }
    # Probably for QPs p*v*v (?)
    quadratic_constraint_variables_associated_to_parameters_cache::Dict{
        MOI.ConstraintIndex,
        Vector{MOI.ScalarAffineTerm{T}},
    }
    # Store the map between quadratic terms add as var * parameter to the resulting shape when implemented in the solver
    # for instance p*p + var -> ScalarAffine(var)
    quadratic_added_cache::OrderedDict{MOI.ConstraintIndex,MOI.ConstraintIndex}
    last_quad_add_added::Int64
    vector_constraint_cache::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.VectorAffineTerm{Float64}},
    }
    # Ditto from the constraint caches but for the objective function
    affine_objective_cache::Vector{MOI.ScalarAffineTerm{T}}
    quadratic_objective_cache_pv::Vector{MOI.ScalarQuadraticTerm{T}}
    quadratic_objective_cache_pp::Vector{MOI.ScalarQuadraticTerm{T}}
    quadratic_objective_cache_pc::Vector{MOI.ScalarAffineTerm{T}}
    quadratic_objective_variables_associated_to_parameters_cache::Vector{
        MOI.ScalarAffineTerm{T},
    }
    multiplicative_parameters::Set{Int64}
    dual_value_of_parameters::Vector{Float64}
    evaluate_duals::Bool
    number_of_parameters_in_model::Int64
    function Optimizer(optimizer::OT; evaluate_duals::Bool = true) where {OT}
        return new{Float64,OT}(
            optimizer,
            Dict{MOI.VariableIndex,Float64}(),
            Dict{MOI.VariableIndex,String}(),
            Dict{MOI.VariableIndex,Float64}(),
            Dict{MOI.VariableIndex,MOI.VariableIndex}(),
            0,
            PARAMETER_INDEX_THRESHOLD,
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.ScalarAffineTerm{Float64}},
            }(),
            Dict{MOI.ConstraintIndex,Vector{MOI.ScalarQuadraticTerm{Float64}}}(),
            Dict{MOI.ConstraintIndex,Vector{MOI.ScalarQuadraticTerm{Float64}}}(),
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.ScalarAffineTerm{Float64}},
            }(),
            Dict{MOI.ConstraintIndex,Vector{MOI.ScalarAffineTerm{Float64}}}(),
            OrderedDict{MOI.ConstraintIndex,MOI.ConstraintIndex}(),
            0,
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.VectorAffineTerm{Float64}},
            }(),
            Vector{MOI.ScalarAffineTerm{Float64}}(),
            Vector{MOI.ScalarQuadraticTerm{Float64}}(),
            Vector{MOI.ScalarQuadraticTerm{Float64}}(),
            Vector{MOI.ScalarAffineTerm{Float64}}(),
            Vector{MOI.ScalarAffineTerm{Float64}}(),
            Set{Int64}(),
            Vector{Float64}(),
            evaluate_duals,
            0,
        )
    end
end

function ParametricOptimizer(args...; kwargs...)
    @warn("ParametricOptimizer is deprecated. Use `POI.Optimizer` instead.")
    return Optimizer(args...; kwargs...)
end

include("utils.jl")
include("duals.jl")
include("update_parameters.jl")

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.optimizer) &&
           isempty(model.parameters) &&
           isempty(model.parameters_name) &&
           isempty(model.variables) &&
           isempty(model.updated_parameters) &&
           isempty(model.variables) &&
           model.last_variable_index_added == 0 &&
           model.last_parameter_index_added == PARAMETER_INDEX_THRESHOLD &&
           isempty(model.affine_constraint_cache) &&
           isempty(model.quadratic_constraint_cache_pv) &&
           isempty(model.quadratic_constraint_cache_pp) &&
           isempty(model.quadratic_constraint_cache_pc) &&
           isempty(
               model.quadratic_constraint_variables_associated_to_parameters_cache,
           ) &&
           isempty(model.quadratic_added_cache) &&
           model.last_quad_add_added == 0 &&
           isempty(model.affine_objective_cache) &&
           isempty(model.quadratic_objective_cache_pv) &&
           isempty(model.quadratic_objective_cache_pp) &&
           isempty(model.quadratic_objective_cache_pc) &&
           isempty(
               model.quadratic_objective_variables_associated_to_parameters_cache,
           ) &&
           isempty(model.dual_value_of_parameters) &&
           model.number_of_parameters_in_model == 0
end

function MOI.supports_constraint(
    model::Optimizer,
    F::Union{
        Type{MOI.SingleVariable},
        Type{MOI.ScalarAffineFunction{T}},
        Type{MOI.VectorOfVariables},
        Type{MOI.VectorAffineFunction{T}},
    },
    S::Type{<:MOI.AbstractSet},
) where {T}
    return MOI.supports_constraint(model.optimizer, F, S)
end

function MOI.supports_constraint(
    model::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{T}},
    S::Type{<:MOI.AbstractSet},
) where {T}
    return MOI.supports_constraint(
        model.optimizer,
        MOI.ScalarAffineFunction{T},
        S,
    )
end

function MOI.supports_constraint(
    model::Optimizer,
    ::Type{MOI.VectorQuadraticFunction{T}},
    S::Type{<:MOI.AbstractSet},
) where {T}
    return MOI.supports_constraint(
        model.optimizer,
        MOI.VectorAffineFunction{T},
        S,
    )
end

function MOI.supports(
    model::Optimizer,
    attr::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction{MOI.SingleVariable},
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}},
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}},
    },
) where {T}
    return MOI.supports(model.optimizer, attr)
end

function MOI.Utilities.supports_default_copy_to(model::Optimizer, bool::Bool)
    return MOI.Utilities.supports_default_copy_to(model.optimizer, bool)
end
function MOI.supports(model::Optimizer, ::MOI.Name)
    return MOI.supports(model.optimizer, MOI.Name())
end
function MOI.get(model::Optimizer, ::MOI.ListOfModelAttributesSet)
    return MOI.get(model.optimizer, MOI.ListOfModelAttributesSet())
end
MOI.get(model::Optimizer, ::MOI.Name) = MOI.get(model.optimizer, MOI.Name())
function MOI.set(model::Optimizer, ::MOI.Name, name::String)
    return MOI.set(model.optimizer, MOI.Name(), name)
end
function MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices)
    return MOI.get(model.optimizer, MOI.ListOfVariableIndices())
end
function MOI.get(model::Optimizer, ::MOI.ListOfVariableAttributesSet)
    return MOI.get(model.optimizer, MOI.ListOfVariableAttributesSet())
end
function MOI.get(model::Optimizer, ::MOI.ListOfConstraintAttributesSet)
    return MOI.get(model.optimizer, MOI.ListOfConstraintAttributesSet())
end
function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{F,S},
    f::F,
) where {F,S}
    MOI.set(model.optimizer, MOI.ConstraintFunction(), c, f)
    return
end
function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{F,S},
    s::S,
) where {F,S}
    MOI.set(model.optimizer, MOI.ConstraintSet(), c, s)
    return
end
function MOI.modify(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
    chg::MOI.ScalarCoefficientChange{Float64},
) where {F,S}
    MOI.modify(model.optimizer, c, chg)
    return
end
function MOI.modify(
    model::Optimizer,
    c::MOI.ObjectiveFunction{F},
    chg::MOI.ScalarCoefficientChange{Float64},
) where {F<:MathOptInterface.AbstractScalarFunction}
    MOI.modify(model.optimizer, c, chg)
    return
end
function MOI.modify(
    model::Optimizer,
    c::MOI.ObjectiveFunction{F},
    chg::MOI.ScalarConstantChange{Float64},
) where {F<:MathOptInterface.AbstractScalarFunction}
    MOI.modify(model.optimizer, c, chg)
    return
end

function MOI.empty!(model::Optimizer{T}) where {T}
    MOI.empty!(model.optimizer)
    empty!(model.parameters)
    empty!(model.parameters_name)
    empty!(model.updated_parameters)
    empty!(model.variables)
    model.last_variable_index_added = 0
    model.last_parameter_index_added = PARAMETER_INDEX_THRESHOLD
    empty!(model.affine_constraint_cache)
    empty!(model.quadratic_constraint_cache_pv)
    empty!(model.quadratic_constraint_cache_pp)
    empty!(model.quadratic_constraint_cache_pc)
    empty!(model.quadratic_constraint_variables_associated_to_parameters_cache)
    empty!(model.quadratic_added_cache)
    model.last_quad_add_added = 0
    empty!(model.vector_constraint_cache)
    empty!(model.affine_objective_cache)
    empty!(model.quadratic_objective_cache_pv)
    empty!(model.quadratic_objective_cache_pp)
    empty!(model.quadratic_objective_cache_pc)
    empty!(model.quadratic_objective_variables_associated_to_parameters_cache)
    empty!(model.dual_value_of_parameters)
    model.number_of_parameters_in_model = 0
    return
end

function MOI.set(
    model::Optimizer,
    attr::MOI.VariableName,
    v::MOI.VariableIndex,
    name::String,
)
    if is_parameter_in_model(model, v)
        model.parameters_name[v] = name
    else
        MOI.set(model.optimizer, attr, v, name)
    end
    return
end

function MOI.get(model::Optimizer, attr::MOI.VariableName, v::MOI.VariableIndex)
    if is_parameter_in_model(model, v)
        return get(model.parameters_name, v, "")
    else
        return MOI.get(model.optimizer, attr, v)
    end
end

function MOI.get(model::Optimizer, tp::Type{MOI.VariableIndex}, attr::String)
    return MOI.get(model.optimizer, tp, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveBound)
    return MOI.get(model.optimizer, attr)
end

function MOI.supports(
    model::Optimizer,
    attr::MOI.VariableName,
    tp::Type{MOI.VariableIndex},
)
    return MOI.supports(model.optimizer, attr, tp)
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ConstraintName,
    c::MOI.ConstraintIndex,
    name::String,
)
    MOI.set(model.optimizer, attr, c, name)
    return
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintName,
    c::MOI.ConstraintIndex,
)
    return MOI.get(model.optimizer, attr, c)
end

function MOI.get(model::Optimizer, ::MOI.NumberOfVariables)
    return length(model.parameters) + length(model.variables)
end

function MOI.supports(
    model::Optimizer,
    attr::MOI.ConstraintName,
    tp::Type{<:MOI.ConstraintIndex},
)
    return MOI.supports(model.optimizer, attr, tp)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    # Check if the index was cached as Affine expression
    if haskey(model.affine_constraint_cache, ci)
        return model.affine_constraint_cache[ci]
        # Check if the index was cached as a Quadratic
    elseif haskey(model.quadratic_added_cache, ci)
        return model.quadratic_added_cache[ci]
    else
        return MOI.get(model.optimizer, attr, ci)
    end
end

function MOI.get(
    model::Optimizer,
    tp::Type{MOI.ConstraintIndex{F,S}},
    attr::String,
) where {F,S}
    return MOI.get(model.optimizer, tp, attr)
end

function MOI.get(model::Optimizer, tp::Type{MOI.ConstraintIndex}, attr::String)
    return MOI.get(model.optimizer, tp, attr)
end

function MOI.is_valid(model::Optimizer, vi::MOI.VariableIndex)
    return MOI.is_valid(model.optimizer, vi)
end

function MOI.supports(model::Optimizer, ::MOI.NumberOfThreads)
    return MOI.supports(model.optimizer, MOI.NumberOfThreads())
end

function MOI.supports(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.supports(model.optimizer, MOI.TimeLimitSec())
end

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)
    MOI.set(model.optimizer, MOI.TimeLimitSec(), value)
    return
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(model.optimizer, MOI.TimeLimitSec())
end

function MOI.get(model::Optimizer, ::MOI.SolveTime)
    return MOI.get(model.optimizer, MOI.SolveTime())
end

function MOI.supports(model::Optimizer, ::MOI.Silent)
    return MOI.supports(model.optimizer, MOI.Silent())
end

function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    MOI.set(model.optimizer, MOI.Silent(), value)
    return
end

MOI.get(model::Optimizer, ::MOI.Silent) = MOI.get(model.optimizer, MOI.Silent())

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return MOI.get(model.optimizer, MOI.RawStatusString())
end

function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{F,S}) where {F,S}
    return length(MOI.get(model, MOI.ListOfConstraintIndices{F,S}()))
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.throw_if_not_valid(model, ci)
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{MOI.SingleVariable,MOI.VectorOfVariables,MOI.VectorAffineFunction},
    S<:MOI.AbstractSet,
}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveSense)
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveFunctionType)
    if !isempty(model.quadratic_constraint_cache_pv) ||
       !isempty(model.quadratic_constraint_cache_pc)
        return MOI.ScalarQuadraticFunction{Float64}
    end
    return MOI.get(model.optimizer, attr)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ObjectiveFunction{F},
) where {
    F<:Union{
        MOI.SingleVariable,
        MOI.ScalarAffineFunction{T},
        MOI.ScalarQuadraticFunction{T},
    },
} where {T}

    # Check if the index was cached as Affine expression
    if !isempty(model.affine_objective_cache)
        return model.affine_objective_cache
        # Check if the index was cached as a Quadratic
    elseif !isempty(model.quadratic_objective_cache_pv)
        return model.quadratic_objective_cache_pv
    elseif !isempty(model.quadratic_objective_cache_pp)
        return model.quadratic_objective_cache_pp
    elseif !isempty(model.quadratic_objective_cache_pc)
        return model.quadratic_objective_cache_pc
    else
        return MOI.get(model.optimizer, attr)
    end
end

function MOI.get(model::Optimizer, attr::MOI.ResultCount)
    return MOI.get(model.optimizer, attr)
end

# In the AbstractBridgeOptimizer, we collect all the possible constraint types and them filter with NumberOfConstraints.
# If NumberOfConstraints is zero then we remove it from the list.
# Here, you can look over keys(quadratic_added_cache) and add the F-S types of all the keys in constraints.
# To implement NumberOfConstraints, you call NumberOfConstraints to the inner optimizer.
# Then you remove the number of constraints of that that in values(quadratic_added_cache)
function MOI.get(model::Optimizer, ::MOI.ListOfConstraints)
    inner_ctrs = MOI.get(model.optimizer, MOI.ListOfConstraints())
    if !has_quadratic_constraint_caches(model)
        return inner_ctrs
    end

    cache_keys = collect(keys(model.quadratic_added_cache))
    constraints = Set{Tuple{DataType,DataType}}()

    for (F, S) in inner_ctrs
        inner_index =
            MOI.get(model.optimizer, MOI.ListOfConstraintIndices{F,S}())
        cache_map_check =
            quadratic_constraint_cache_map_check(mode, inner_index)
        push!(constraints, typeof.(cache_keys[cache_map])...)
        # If not all the constraints are chached then also push the original type
        if !all(cache_map_check)
            push!(constraints, (F, S))
        end
    end

    return collect(constraints)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ListOfConstraintIndices{F,S},
) where {S,F<:Union{MOI.VectorOfVariables,MOI.SingleVariable}}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{F,S},
) where {
    S<:MOI.AbstractSet,
    F<:Union{MOI.ScalarAffineFunction{T},MOI.VectorAffineFunction{T}},
} where {T}
    inner_index = MOI.get(model.optimizer, MOI.ListOfConstraintIndices{F,S}())
    if !has_quadratic_constraint_caches(model)
        return inner_index
    end

    cache_map_check = quadratic_constraint_cache_map_check(mode, inner_index)
    return inner_index[cache_map_check]
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{F,S},
) where {S<:MOI.AbstractSet,F<:MOI.ScalarQuadraticFunction{T}} where {T}
    inner_index = MOI.ConstraintIndex{F,S}[]
    if MOI.supports_constraint(model.optimizer, F, S)
        inner_index =
            MOI.get(model.optimizer, MOI.ListOfConstraintIndices{F,S}())
        if !has_quadratic_constraint_caches(model)
            return inner_index
        end
    end

    quadratic_caches = [
        :quadratic_constraint_cache_pc,
        :quadratic_constraint_cache_pp,
        :quadratic_constraint_cache_pv,
        # JD: Check if this applies here
        # :quadratic_constraint_variables_associated_to_parameters_cache
    ]

    for field in quadratic_caches
        cache = getfield(model, field)
        push!(inner_index, keys(cache)...)
    end
    return inner_index
end

function MOI.supports_add_constrained_variable(::Optimizer, ::Type{Parameter})
    return true
end

function MOI.supports_add_constrained_variables(
    model::Optimizer,
    ::Type{MOI.Reals},
)
    return MOI.supports_add_constrained_variables(model.optimizer, MOI.Reals)
end

function MOI.add_variable(model::Optimizer)
    next_variable_index!(model)
    v_p = MOI.VariableIndex(model.last_variable_index_added)
    v = MOI.add_variable(model.optimizer)
    model.variables[v_p] = v
    return v_p
end

function MOI.add_constrained_variable(model::Optimizer, set::Parameter)
    next_parameter_index!(model)
    p = MOI.VariableIndex(model.last_parameter_index_added)
    model.parameters[p] = set.val
    cp = MOI.ConstraintIndex{MOI.SingleVariable,Parameter}(
        model.last_parameter_index_added,
    )
    update_number_of_parameters!(model)
    return p, cp
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.SingleVariable,
    set::MOI.AbstractScalarSet,
)
    if is_parameter_in_model(model, f.variable)
        error("Cannot constrain a parameter")
    elseif !is_variable_in_model(model, f.variable)
        error("Variable not in the model")
    end
    return MOI.add_constraint(model.optimizer, f, set)
end

function add_constraint_with_parameters_on_function(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
    set::MOI.AbstractScalarSet,
) where {T}
    vars, params, param_constant =
        separate_possible_terms_and_calculate_parameter_constant(model, f.terms)
    ci = MOI.Utilities.normalize_and_add_constraint(
        model.optimizer,
        MOI.ScalarAffineFunction(vars, f.constant + param_constant),
        set,
    )
    model.affine_constraint_cache[ci] = params
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
    set::MOI.AbstractScalarSet,
) where {T}
    if !function_has_parameters(model, f)
        return MOI.add_constraint(model.optimizer, f, set)
    else
        return add_constraint_with_parameters_on_function(model, f, set)
    end
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    v::MOI.VariableIndex,
)
    if is_parameter_in_model(model, v)
        return model.parameters[v]
    elseif is_variable_in_model(model, v)
        return MOI.get(model.optimizer, attr, model.variables[v])
    else
        error("Variable not in the model")
    end
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintSet,
    cp::MOI.ConstraintIndex{MOI.SingleVariable,Parameter},
    set::Parameter,
)
    p = MOI.VariableIndex(cp.value)
    if !is_parameter_in_model(model, p)
        error("Parameter not in the model")
    end
    return model.updated_parameters[p] = set.val
end

struct ParameterValue <: MOI.AbstractVariableAttribute end

"""
    MOI.set(
        model::Optimizer,
        ::ParameterValue,
        x::MOI.VariableIndex,
        value::Float64,
    )

Sets the parameter to a given value, using its `MOI.VariableIndex` as reference.

# Example

```julia-repl
julia> MOI.set(model, ParameterValue(), w, 2.0)
2.0
```
"""
function MOI.set(
    model::Optimizer,
    ::ParameterValue,
    vi::MOI.VariableIndex,
    val::Float64,
)
    if !is_parameter_in_model(model, vi)
        error("Parameter not in the model")
    end
    return model.updated_parameters[vi] = val
end

function MOI.set(
    model::Optimizer,
    ::ParameterValue,
    vi::MOI.VariableIndex,
    val::Real,
)
    return MOI.set(model, ParameterValue(), vi, convert(Float64, val))
end

function empty_objective_function_caches!(model::Optimizer)
    empty!(model.affine_objective_cache)
    empty!(model.quadratic_objective_cache_pv)
    empty!(model.quadratic_objective_cache_pp)
    empty!(model.quadratic_objective_cache_pc)
    return
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ObjectiveFunction,
    f::MOI.ScalarAffineFunction{T},
) where {T}
    # clear previously defined objetive function cache
    empty_objective_function_caches!(model)
    if !function_has_parameters(model, f)
        MOI.set(model.optimizer, attr, f)
    else
        vars, params, param_constant =
            separate_possible_terms_and_calculate_parameter_constant(
                model,
                f.terms,
            )
        MOI.set(
            model.optimizer,
            attr,
            MOI.ScalarAffineFunction(vars, f.constant + param_constant),
        )
        model.affine_objective_cache = params
    end
    return
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ObjectiveFunction,
    v::MOI.SingleVariable,
)
    if haskey(model.parameters, v)
        error("Cannot use a parameter as objective function alone")
    elseif !haskey(model.variables, v)
        error("Variable not in the model")
    end
    return MOI.set(
        model.optimizer,
        attr,
        MOI.SingleVariable(model.variables[v.variable]),
    )
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    MOI.set(model.optimizer, attr, sense)
    return
end

function MOI.get(
    model::Optimizer,
    attr::T,
) where {
    T<:Union{
        MOI.TerminationStatus,
        MOI.ObjectiveValue,
        MOI.DualObjectiveValue,
        MOI.PrimalStatus,
        MOI.DualStatus,
    },
}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(
    model::Optimizer,
    attr::T,
    c::MOI.ConstraintIndex,
) where {T<:Union{MOI.ConstraintPrimal,MOI.ConstraintDual}}
    return MOI.get(model.optimizer, attr, c)
end

function MOI.set(model::Optimizer, attr::MOI.RawParameter, val::Any)
    MOI.set(model.optimizer, attr, val)
    return
end

# TODO(odow): remove this.
function MOI.set(model::Optimizer, attr::String, val::Any)
    return MOI.set(model.optimizer, MOI.RawParameter(attr), val)
end

function MOI.get(model::Optimizer, ::MOI.SolverName)
    name = MOI.get(model.optimizer, MOI.SolverName())
    return "Parametric Optimizer with $(name) attached"
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VectorOfVariables,
    set::MOI.AbstractVectorSet,
) where {T}
    if function_has_parameters(model, f)
        error("VectorOfVariables does not allow parameters")
    end
    return MOI.add_constraint(model.optimizer, f, set)
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VectorAffineFunction{T},
    set::MOI.AbstractVectorSet,
) where {T}
    if !function_has_parameters(model, f)
        return MOI.add_constraint(model.optimizer, f, set)
    else
        return add_constraint_with_parameters_on_function(model, f, set)
    end
end

function add_constraint_with_parameters_on_function(
    model::Optimizer,
    f::MOI.VectorAffineFunction{T},
    set::MOI.AbstractVectorSet,
) where {T}
    vars, params, param_constants =
        separate_possible_terms_and_calculate_parameter_constant(model, f, set)
    ci = MOI.add_constraint(
        model.optimizer,
        MOI.VectorAffineFunction(vars, f.constants + param_constants),
        set,
    )
    model.vector_constraint_cache[ci] = params
    return ci
end

function add_constraint_with_parameters_on_function(
    model::Optimizer,
    f::MOI.ScalarQuadraticFunction{T},
    set::S,
) where {T,S<:MOI.AbstractScalarSet}
    if function_quadratic_terms_has_parameters(model, f.quadratic_terms)
        (
            quad_vars,
            quad_aff_vars,
            quad_params,
            aff_terms,
            variables_associated_to_parameters,
            quad_param_constant,
        ) = separate_possible_terms_and_calculate_parameter_constant(
            model,
            f.quadratic_terms,
        )
    else
        quad_vars = f.quadratic_terms
        quad_aff_vars = MOI.ScalarQuadraticTerm{T}[]
        quad_params = MOI.ScalarQuadraticTerm{T}[]
        aff_terms = MOI.ScalarAffineTerm{T}[]
        quad_param_constant = zero(T)
        variables_associated_to_parameters = MOI.VariableIndex[]
    end
    if function_affine_terms_has_parameters(model, f.affine_terms)
        (
            aff_vars,
            aff_params,
            terms_with_variables_associated_to_parameters,
            aff_param_constant,
        ) = separate_possible_terms_and_calculate_parameter_constant(
            model,
            f.affine_terms,
            variables_associated_to_parameters,
        )
    else
        aff_vars = f.affine_terms
        aff_params = MOI.ScalarAffineTerm{T}[]
        terms_with_variables_associated_to_parameters =
            MOI.ScalarAffineTerm{T}[]
        aff_param_constant = zero(T)
    end
    aff_terms = vcat(aff_terms, aff_vars)
    const_term = f.constant + aff_param_constant + quad_param_constant
    quad_terms = quad_vars
    f_quad = if !isempty(quad_vars)
        MOI.ScalarQuadraticFunction(aff_terms, quad_terms, const_term)
    else
        MOI.ScalarAffineFunction(aff_terms, const_term)
    end
    model.last_quad_add_added += 1
    ci =
        MOI.Utilities.normalize_and_add_constraint(model.optimizer, f_quad, set)
    # This part is used to remember that ci came from a quadratic function
    # It is particularly useful because sometimes the constraint mutates
    new_ci = MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{T},S}(
        model.last_quad_add_added,
    )
    model.quadratic_added_cache[new_ci] = ci
    fill_quadratic_constraint_caches!(
        model,
        new_ci,
        quad_aff_vars,
        quad_params,
        aff_params,
        terms_with_variables_associated_to_parameters,
    )
    return new_ci
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.ScalarQuadraticFunction{T},
    set::MOI.AbstractScalarSet,
) where {T}
    if !function_has_parameters(model, f)
        return MOI.add_constraint(model.optimizer, f, set)
    else
        return add_constraint_with_parameters_on_function(model, f, set)
    end
end

function MOI.delete(model::Optimizer, v::MOI.VariableIndex)
    pop!(model.variables, v, 0) # 0 is a default value if v is not found
    MOI.delete(model.optimizer, v)
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
) where {F<:MOI.ScalarAffineFunction,S<:MOI.AbstractSet}
    if haskey(model.affine_constraint_cache, c)
        delete!(model.affine_constraint_cache, c)
    end
    MOI.delete(model.optimizer, c)
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{MOI.SingleVariable,MOI.VectorOfVariables,MOI.VectorAffineFunction},
    S<:MOI.AbstractSet,
}
    MOI.delete(model.optimizer, c)
    return
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{MOI.SingleVariable,MOI.VectorOfVariables,MOI.VectorAffineFunction},
    S<:MOI.AbstractSet,
}
    return MOI.is_valid(model.optimizer, c)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
) where {F<:MOI.ScalarAffineFunction,S<:MOI.AbstractSet}
    if haskey(model.affine_constraint_cache, c)
        # TODO: how to check the cached constraint is valid (?)
        return function_has_parameters(
            model,
            model.affine_constraint_cache[c],
        ) && MOI.is_valid(model.optimizer, c)
    end
    return MOI.is_valid(model.optimizer, c)
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ObjectiveFunction,
    f::MOI.ScalarQuadraticFunction{T},
) where {T}
    # clear previously defined objetive function cache
    empty_objective_function_caches!(model)
    if !function_has_parameters(model, f)
        MOI.set(model.optimizer, attr, f)
        return
    end
    if function_quadratic_terms_has_parameters(model, f.quadratic_terms)
        (
            quad_vars,
            quad_aff_vars,
            quad_params,
            aff_terms,
            variables_associated_to_parameters,
            quad_param_constant,
        ) = separate_possible_terms_and_calculate_parameter_constant(
            model,
            f.quadratic_terms,
        )
    else
        quad_vars = f.quadratic_terms
        quad_aff_vars = MOI.ScalarQuadraticTerm{T}[]
        quad_params = MOI.ScalarQuadraticTerm{T}[]
        aff_terms = MOI.ScalarAffineTerm{T}[]
        quad_param_constant = zero(T)
        variables_associated_to_parameters = MOI.VariableIndex[]
    end

    if function_affine_terms_has_parameters(model, f.affine_terms)
        (
            aff_vars,
            aff_params,
            terms_with_variables_associated_to_parameters,
            aff_param_constant,
        ) = separate_possible_terms_and_calculate_parameter_constant(
            model,
            f.affine_terms,
            variables_associated_to_parameters,
        )
    else
        aff_vars = f.affine_terms
        aff_params = MOI.ScalarAffineTerm{T}[]
        terms_with_variables_associated_to_parameters =
            MOI.ScalarAffineTerm{T}[]
        aff_param_constant = zero(T)
    end

    aff_terms = vcat(aff_terms, aff_vars)
    const_term = f.constant + aff_param_constant + quad_param_constant
    quad_terms = quad_vars

    if !isempty(quad_vars)
        MOI.set(
            model.optimizer,
            attr,
            MOI.ScalarQuadraticFunction(aff_terms, quad_terms, const_term),
        )
    else
        MOI.set(
            model.optimizer,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
            MOI.ScalarAffineFunction(aff_terms, const_term),
        )
    end

    if !isempty(quad_terms)
        f_quad = MOI.ScalarQuadraticFunction(aff_terms, quad_terms, const_term)

        MOI.set(model.optimizer, attr, f_quad)
    else
        f_quad = MOI.ScalarAffineFunction(aff_terms, const_term)

        MOI.set(
            model.optimizer,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
            f_quad,
        )
    end

    model.quadratic_objective_cache_pv = quad_aff_vars
    model.quadratic_objective_cache_pp = quad_params
    model.quadratic_objective_cache_pc = aff_params
    model.quadratic_objective_variables_associated_to_parameters_cache =
        terms_with_variables_associated_to_parameters

    return
end

function MOI.optimize!(model::Optimizer)
    if !isempty(model.updated_parameters)
        update_parameters!(model)
    end
    MOI.optimize!(model.optimizer)
    if model.evaluate_duals
        calculate_dual_of_parameters(model)
    end
    return
end

end # module


text:
abstract type AbstractCircularPolygon <: AbstractClosedPath end
abstract type AbstractPolygon <: AbstractCircularPolygon end

function show(io::IO,P::AbstractCircularPolygon)
	print(IOContext(io,:compact=>true),typeof(P)," with ",length(P)," sides")
end
function show(io::IO,::MIME"text/plain",P::AbstractCircularPolygon)
	print(io,typeof(P)," with ",length(P)," sides")
end

# Other methods
sides(p::AbstractCircularPolygon) = curves(p)
side(p::AbstractCircularPolygon,args...) = curve(p,args...)

function reverse(p::AbstractCircularPolygon)
	@assert sum(isinf.(vertices(p))) < 2 "Reversal is not a propoer polygon"
	typeof(p)(reverse(reverse.(sides(p))))
end

function winding(p::AbstractCircularPolygon,z::Number)
	if isinf(z)
		!isfinite(p) && @warn "Winding number is ill-defined for a boundary point"
		return 0
	end

	# truncate an unbounded path
	if !isfinite(p)
		C = enclosing_circle(p)
		while !isinside(z,C)
			C = Circle(C.center,2*C.radius)
		end
		return winding(truncate(p,C),z)
	end

	w = 0
	v = vertex(p,1)
	for s in p
		vnew = point(s,1)
		if s isa Segment
			w += angle((vnew-z)/(v-z))
		else  #if s isa Arc
			# move the branch cut to avoid the arc
			if abs(z-s.circle.center) < s.circle.radius
				# can use ray from the center to a point not on s
				u = point(s.circle,s.start-(1-s.delta)/2) # not on s
				w += angle((vnew-z)/(z-u)) - angle((v-z)/(z-u))
			else
				# can put center on the positive real axis
				u = s.circle.center
				w += angle((vnew-z)/(u-z)) - angle((v-z)/(u-z))
			end
		end
		v = vnew
	end
	return round(Int,w/(2π))
end

#
# CircularPolygon
#
"""
	(type) CircularPolygon
Type for closed paths consisting entirely of arcs, segments, and rays.
"""
struct CircularPolygon <: AbstractCircularPolygon
	path
	function CircularPolygon(p::AbstractClosedPath)
		# Continuity and closure have been checked to make a closed path
		valid = isa.(curves(p),Union{Arc,Segment,Ray})
		@assert all(valid) "All sides must be an Arc, Segment, or Ray"
		new(p)
	end
end

# Constructors
"""
	CircularPolygon(p::AbstractPath; tol=<default>)
	CircularPolygon(p::AbstractVector; tol=<default>)
Construct a circular polygon from a (possibly closed) path, or from a vector of curves. The `tol` parameter is a tolerance used when checking continuity and closedness of the path.
"""
CircularPolygon(p::AbstractPath;kw...) = CircularPolygon(ClosedPath(p;kw...))
function CircularPolygon(p::AbstractVector{T};kw...) where T<:AbstractCurve
	CircularPolygon(ClosedPath(p;kw...))
end

# Required methods
curves(p::CircularPolygon) = curves(p.path)
curve(p::CircularPolygon,k::Integer) = curve(p.path,k)
arclength(p::CircularPolygon) = arclength(p.path)
(p::CircularPolygon)(t) = point(p.path,t)

inv(p::CircularPolygon) = CircularPolygon([inv(c) for c in curves(p)])

# Other methods
"""
	ispositive(p::CircularPolygon)
Determine whether the circular polygon is positively oriented (i.e., circulates counterclockwise around the points it encloses).
"""
ispositive(p::CircularPolygon) = winding(0,1/p) < 0

#
# Polygon
#

# Type
"""
	(type) Polygon
Type for closed paths consisting entirely of segments and rays.
"""
struct Polygon <: AbstractPolygon
	path
	function Polygon(p::AbstractClosedPath) where T<:AbstractCurve
		# Assumes continuity and closure have been checked previously
		valid = isa.(curves(p),Union{Segment,Ray})
		@assert all(valid) "All sides must be a Segment or Ray"
		new(p)
	end
end

# Constructors
"""
	Polygon(p::AbstractPath; tol=<default>)
	Polygon(p::AbstractVector{T<:AbstractCurve}; tol=<default>)
Construct a polygon from a (possibly closed) path, or from a vector of curves. The `tol` parameter is a tolerance used when checking continuity and closedness of the path.
"""
Polygon(p::AbstractPath;kw...) = Polygon(ClosedPath(p;kw...))
function Polygon(p::AbstractVector{T};kw...) where T<:AbstractCurve
	Polygon(ClosedPath(p;kw...))
end

"""
	Polygon(v::AbstractVector)
Construct a polygon from a vector of its vertices. Each element of `v` should be either a finite vertex, or a tuple of two angles that indicate the angles of two rays incident to an infinite vertex: one "to" infinity, and a second "from" infinity.
"""
function Polygon(v::AbstractVector)
	n = length(v)
	p = Vector{Union{Segment,Ray}}(undef,n)
	for j = 1:n
		vthis = v[j]
		vnext = v[mod(j,n)+1]
		if isa(vthis,Tuple)
			if isa(vnext,Tuple)
				@error("Cannot have consecutive infinite vertices")
				return nothing
			else
				p[j] = Ray(vnext,vthis[2],true)
			end
		else
			if isa(vnext,Tuple)
				p[j] = Ray(vthis,vnext[1])
			else
				p[j] = Segment(vthis,vnext)
			end
		end
	end
	return Polygon(ClosedPath(p))
end

# Required methods
curves(p::Polygon) = curves(p.path)
arclength(p::Polygon) = arclength(p.path)
(p::Polygon)(t) = point(p.path,t)

inv(p::Polygon) = CircularPolygon([inv(c) for c in curves(p)])

# Display methods
function show(io::IO,::MIME"text/plain",P::Polygon)
	print(io,"Polygon with ",length(P)," vertices:")
	for (v,a) in zip(vertices(P),angles(P))
		print(io,"\n   ")
		show(io,MIME("text/plain"),v)
		print(io,", interior angle ",a/pi,"⋅π")
	end
end

# Other methods
"""
	truncate(P::Union{CircularPolygon,Polygon})
Apply `truncate` to `P` using a circle that is centered at the centroid of its finite vertices, and a radius twice the maximum from the centroid to the finite vertices.
"""
function truncate(p::Union{CircularPolygon,Polygon})
	isfinite(p) && return p   # nothing to do
	return truncate(p,enclosing_circle(p,4))
end

"""
	truncate(P::Union{CircularPolygon,Polygon},C::Circle)
Compute a trucated form of the polygon by replacing each pair of rays incident at infinity with two segments connected by an arc along the given circle. This is *not* a true clipping of the polygon, as finite sides are not altered. The result is either a CircularPolygon or the original `P`.
"""
function truncate(p::Union{CircularPolygon,Polygon},c::Circle)
	n = length(p)
	s,v = sides(p),vertices(p)
	snew = Vector{Any}(undef,n)
	z_pre = NaN
	if isa(s[1],Ray) && isinf(point(s[1],0))
		# recognize that first side is actually the return from an infinite vertex
		z_pre = intersect(c,s[n])[1]
	end
	for k = 1:n
		if !isa(s[k],Ray)
			snew[k] = s[k]
		else
			# first of a pair?
			if isnan(z_pre)
				z_pre = intersect(c,s[k])[1]
				snew[k] = Segment(s[k].base,z_pre)
			else
				z_post = intersect(c,s[k])[1]
				t0 = arg(c,z_pre)
				delta = mod( angle( (z_post-c.center)/(z_pre-c.center) ) / (2π), 1)
				snew[k] = [Arc(c,t0,delta),Segment(z_post,s[k].base)]
				z_pre = NaN
			end
		end
	end
	return CircularPolygon(vcat(snew...))
end

"""
	angles(P::Polygon)
Compute a vector of interior angles at the vertices of the polygon `P`. At a finite vertex these lie in (0,2π]; at an infinite vertex, the angle is in [-2π,0].
"""
function angles(p::Polygon)
	# computes a turn angle in (-pi,pi]  (neg = left turn)
	turn(s1,s2) = π - mod(angle(s2/s1)+π,2π)
	s = unittangent.(p)
	n = length(p)
	v = vertices(p)
	θ = similar(real(s))
	for k = 1:n
		kprev = mod(k-2,n)+1
		θ[k] = π + turn(s[kprev],s[k])
		if isinf(v[k])
			θ[k] -= 2π
			if θ[k]==0
				# need a finite perturbation to distinguish 0,-2
				R = maximum(abs.(filter(isfinite,v)))
				C = Circle(0,100*R)
				zprev = intersect(p[kprev],C)
				znext = intersect(p[k],C)
				if angle(znext[1]/zprev[1]) < 0
					θ[k] -= 2π
				end
			end
		end
	end
	# # correct for possible clockwise orientation
	# sum(θ) > 0 ? θ : -θ
	return θ
end

"""
	ispositive(p::Polygon)
Determine whether the polygon is positively oriented (i.e., circulates counterclockwise around the points it encloses).
"""
ispositive(p::Polygon) = sum(angles(p)/pi .- 1) < 0

# Special polygon constructors

"""
	rectangle(xlim,ylim)
Construct the rectangle defined by `xlim[1]`` < x < `xlim[2]`, `ylim[1]`` < y < `ylim[2]`, where z=complex(x,y).
"""
function rectangle(xlim::AbstractVector,ylim::AbstractVector)
	x = [xlim[1],xlim[2],xlim[2],xlim[1],xlim[1]]
	y = [ylim[1],ylim[1],ylim[2],ylim[2],ylim[1]]
	Polygon( [Segment(complex(x[k],y[k]),complex(x[k+1],y[k+1])) for k in 1:4] )
end

"""
	rectangle(z1,z2)
Construct the rectangle whose opposing corners are the given complex values.
"""
rectangle(z1::AnyComplex,z2::AnyComplex) = rectangle(sort(real([z1,z2])),sort(imag([z1,z2])))
rectangle(z1::Number,z2::Number) = rectangle(promote(complex(float(z1)),complex(float(z2)))...)

"""
	n_gon(n)
Construct a regular n-gon with vertices on the unit circle.
"""
function n_gon(n::Integer)
	@assert n > 2 "Must have at least three vertices"
	Polygon( exp.(2im*pi*(0:n-1)/n) )
end