FittingInterface.cpp

#include <iostream>
#include <string>
#include <utility>
#include <iomanip>
#include <vector>
#include <variant>    // for std::variant
#include "ManySimulations.cpp"
#include "cellLibrary/ModelA549.h"


class FittingInterface {
public:
    // PROPERTIES
    double maxDose;
    std::vector<ProtractedDamageInput> scheduledProtractedDamageInputCell;
    std::string checkManySimulationStopMode = "fast";

    // CONSTRUCTOR: load the SDD once
    FittingInterface(const std::string& sddPath, double maxDose = 2)
            : manySdd(sddPath, maxDose), maxDose(maxDose)
    {
//        std::cerr << "[FittingInterface] SDD loaded from: " << sddPath << std::endl;
    }

    // RUN FUNCTIONS
    // Run simulation with given RepairMechanism parameters
    // Return pair of vectors: (times, meanNumberDamages)
    std::pair<std::vector<double>, std::vector<double>> runSimulationFoci(
            double detectRate,
            double finishRate,
            double correctProb
    ) {
        // Create Cell
        Cell cell(manySdd, 0, 2, 0, -1, false); // MUST BE FALSE. MANY SIMS REQUIRES IT
        cell.repairModel = RepairMechanism(detectRate,finishRate, correctProb);
        // Run ManySimulations
        int nSims = 100;
        ManySimulations manySims(cell, nSims);
        manySims.manyRun(60*25, 15, false);
        // Compute average number of damage sites over time
        auto [timePoints, mean, stddev, stderr] =
                computeStatsOverTime(
                        manySims.manyArraySimulationTime,
                        manySims.manyArrayNumberDetectedDamageSitesPerTime, true
                );
        return { timePoints, mean };
    }
    // Return survival:
    double runSimulationSurvival(
            int thresholdForDeath,
            int thresholdCycleChange,
            double weightLethalityMisrepair
    ) {
        // Create Cell
        Cell cell(manySdd, 0, 2, 0, -1, false); // MUST BE FALSE. MANY SIMS REQUIRES IT
        cell.repairModel = RepairMechanism(0.2, 0.0022, -0.30756086282);
        // Variables to fit
        cell.thresholdForDeath = thresholdForDeath;
        cell.thresholdCycleChange = thresholdCycleChange;
        cell.weightLethalityMisrepair = weightLethalityMisrepair;
        // Run ManySimulations
        int nSims = 10000;
        ManySimulations manySims(cell, nSims);
        manySims.manyRun(60*25, 15, false);
        // Compute average number of damage sites over time
        manySims.calculateAverageSurvivalAndError();
        return manySims.averageSurvival;
    }

    // Return survival:
    double runSimulationSurvivalAllParameters(
            double detectRate,
            double finishRate,
            double correctProb,
            int thresholdForDeath,
            int thresholdCycleChange,
            double weightLethalityMisrepair
    ) {
        // Create Cell
        Cell cell(manySdd, 0, 2, 0, -1, false); // MUST BE FALSE. MANY SIMS REQUIRES IT
        cell.repairModel = RepairMechanism(detectRate, finishRate, correctProb);
        cell.thresholdForDeath = thresholdForDeath;
        cell.thresholdCycleChange = thresholdCycleChange;
        cell.weightLethalityMisrepair = weightLethalityMisrepair;
        // Run ManySimulations
        int nSims = 100;
        ManySimulations manySims(cell, nSims);
        manySims.manyRun(60*25, 15, false);
        // Compute average number of damage sites over time
        manySims.calculateAverageSurvivalAndError();
        return manySims.averageSurvival;
    }

    // Generalized simulation runner
    std::variant<
    std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>, // for foci (time, mean, stderr)
    std::pair<double, double>                                                  // for survival (mean, stderr)
> runSimulationGeneral(
            const std::string& mode, // "foci" or "survival"
            int nSims,
            double totalTime,
            double timeStep,
            double detectRate,
            double finishRate,
            double correctProb,
            int thresholdForDeath,
            int thresholdCycleChange,
            double weightLethalityMisrepair,
            int thresholdForNecrosis = 1000,
            double constMisrepairPerNumberOfDamage = -1,
            double newDose = -1
    ) {
        // boolean for daughters
        bool simulateDaughters = false;

        // Create Cell
        // ModelA549 myModel; // can be used
        Cell cell(manySdd, 0, 2, 0, -1,
            false, newDose);
        // Add protracted damages if needed
        if (!scheduledProtractedDamageInputCell.empty()) {
            auto [times, sddPaths, doses] = ProtractedDamageInput::toVectors(scheduledProtractedDamageInputCell);
            auto protractedEvents = ProtractedDamage::loadProtractedDamageSites(times, sddPaths, doses, cell);
            cell.protractedDamageEvents = std::move(protractedEvents);
            // The simulation now lasts totalTime (usually 24-25h) plus the time when you last added the damage
            double maxTime = *std::max_element(times.begin(), times.end());
            // std::cerr << totalTime << std::endl; //debug
            totalTime += maxTime;
            // std::cerr << totalTime << std::endl; //debug
            // ProtractedDamage::printSummary(cell.protractedDamageEvents); //debug
            // ProtractedDamage::printEventSchedule(cell.protractedDamageEvents); //debug
            simulateDaughters = true; //in the future
        }
        cell.repairModel = RepairMechanism(detectRate, finishRate, correctProb);
        cell.thresholdForDeath = thresholdForDeath;
        cell.thresholdCycleChange = thresholdCycleChange;
        cell.weightLethalityMisrepair = weightLethalityMisrepair;
        cell.thresholdForNecrosis = thresholdForNecrosis;
        if (constMisrepairPerNumberOfDamage > 0) {
            cell.boolVariableMisrepair = true;
            cell.repairModel.misrepairMechanism.constMisrepairPerNumberOfDamage = constMisrepairPerNumberOfDamage;
        }

        // Run simulations
        ManySimulations manySims(cell, nSims);
        manySims.checkManySimulationStopMode = checkManySimulationStopMode;
        // std::cerr << "Cell dose: " << cell.doseInCell << std::endl; //debug
        // std::cerr << cell.repairModel.repairDetectionRate << " " << cell.repairModel.repairSpeedRate
        // << " " << cell.repairModel.successRepairConstant << " " << cell.thresholdForDeath
        // << " " << cell.thresholdCycleChange << " " << cell.weightLethalityMisrepair << std::endl; //debug
        manySims.manyRun(totalTime, timeStep, simulateDaughters);

        if (mode == "foci") {
            auto [timePoints, mean, stddev, stderr] =
                computeStatsOverTime(
                    manySims.manyArraySimulationTime,
                    manySims.manyArrayNumberDetectedDamageSitesPerTime, true
                );
            return std::make_tuple(timePoints, mean, stderr);
        }
        else if (mode == "survival") {
            manySims.calculateAverageSurvivalAndError();
            return std::make_pair(manySims.averageSurvival, manySims.stderrSurvival);
        }
        else {
            throw std::invalid_argument("Invalid mode. Use 'foci' or 'survival'.");
        }
    }

private:
    ManySddReader manySdd;
};