---

main.cpp

Go to the documentation of this file.

// CCancer.cpp : Defines the entry point for the console application.
//

#ifdef _WINDOWS
#include "StdAfx.h"
#include <crtdbg.h> // For debugging functions windows only.
#endif

// debugging only
//#define _CRT_DEBUG_MAP_ALLOC

#include "QuorumMaster.h"
#include "InfrastructureMaster.h"

/*#include <fstream.h>
#include <stdlib.h>
#include <iostream.h>
#include <string.h>
//#include "PC12EGFNGFMinimizableDirector.h"
//#include "OneGeneInhibitoryNetworkMinimizableDirector.h"
//#include "AlternateOneGeneInhibitoryNetworkMinimizableDirector.h"
#include "ZeroTAnnealMinimizer.h"
#include "PositiveDefiniteLevenbergMarquardtMinimizer.h"
#include "ParameterReader.h"
#include "GnuPlotterSimulatedAnnealingObserver.h"
#include "LogParameterFilter.h"
#include "SquareRootParameterFilter.h"
#include "IdentityParameterFilter.h"
#include "MixedReactionNetwork.h"
//#include "OneGeneInhibitoryAveragingDirector.h"
//#include "HessianMonteCarloQuench.h"
#include "BasicQuenchMinimizer.h"
#include "RobustLevenbergMarquardtMinimizer.h"
#include "ImprovedLevenbergMarquardtMinimizer.h"
#include "ScaleInvariantMarquardtMinimizer.h"
#include "ConjugateGradientMinimizer.h"
//#include "PC12EGFNGFRunDirector.h"
//#include "PC12AlternateMinimizableDirector.h"
#include "ParameterRandomizer.h"
#include "ExampleNetworkMinimizableDirector.h"
#include "StochasticSensitivityAnalysis.h"
#include "GnuPlotterCostObserver.h"
#include "GnuPlotterParameterScatterPlot.h"
#include "GnuPlotterResidualScatterPlot.h"
#include "CommaStrategyOne.h"
#include "ClonalCommaES.h"
//#include "OneGeneEnsembleRunDirector.h"
//#include "OneGeneEnsembleDirector.h"
//#include "PC12EGFNGFEnsembleDirector.h"
//#include "PC12AlternateEnsembleDirector.h"
//#include "CoupledOneGeneNetworksMinimizableDirector.h"
//#include "CoupledPC12NetworksMinimizableDirector.h"
//#include "EGFRTraffickingCoupledMinimizableDirector.h"
#include "LowDimensionalCostFunctionMapper.h"
#include "QuorumSensingNetworkMinimizableDirector.h"
#include "QuorumSensingRunDirector.h"
//#include "OneGeneInhibitoryNetworkRunDirector.h"
#include "LeastSquaresObserver.h"
#include "ConjugateGradientObserver.h"
//#include "OneGeneStochasticNetworkMinimizableDirector.h"
//#include "LacNetwork.h"
//#include "LacRunDirector.h"
#include "MichaelisMentenHeterodimerizationReaction.h"
//#include "CombinatorialNetwork.h"
//#include "CombinatorialRunDirector.h"
//#include "CombinatorialNetworkMinimizableDirector.h"
#include "LevenbergMarquardtMinimizer.h"
#include "QuorumSensingEnsembleDirector.h"
#include <fenv.h>
*/

main(int argc, char* argv[])
{
/*      int fe=feenableexcept(FE_ALL_EXCEPT);
        if (fe == -1) {
                cout << "feenableexcept unsuccessful\n";
        }
*/      int i;
        std::string defaultkey = "default";
        char *parameterfile;
        char *method;
        int smethod; //to be used in switch
        int plotornot;
        if (argc == 1) {
                parameterfile = "goodtemp6.par";
                method = "c";
                smethod = method[0];
                plotornot = 1;
        }
        else {
                parameterfile = argv[1];
//since switch can only use single characters the code for which method to 
// use is:
// c = conjugate gradient
// l = regular levenberg marquardt
// s = scale invariant levenberg marquardt
// b = basic quench (steepest descent with step size 1.0)
// p = plot (steepest descent with step size 0.0)
// h = stochastic sensitivity analysis (h for Hessian)
// e = ensembledirector (adds error bars to plots)
                if (strlen(argv[2]) == 0){
                        method = "p";
                }
                else {
                        method = argv[2];
                }       
                smethod = method[0];
//if you want to watch the plots, the third argument should be y
                if (strlen(argv[3]) == 0){
                        plotornot = 0;
                } // end if no argv3
                else {
                        if (argv[3][0] == 'y'){
                                plotornot = 1;
                        } // end if argv3 = y
                        else {
                                plotornot = 0;
                        }//end else argv3 = y
                } // end else no argv3
        } // end else strlen

        if (smethod == 'p') {
                plotornot = 1;
        }

//      cout << "plotornot " << plotornot << endl;

        // MICROSOFT DEBUGGING STUFF
        //_crtBreakAlloc = 55;  // Breaks at memory allocation failure no. ?
        // global variable
        //_CrtSetBreakAlloc(1129148);

//      int tmpFlag = _CRTDBG_REPORT_FLAG;
//      tmpFlag |= _CRTDBG_DELAY_FREE_MEM_DF;
//      tmpFlag |= _CRTDBG_ALLOC_MEM_DF; // use any checking at all.
//      tmpFlag |= _CRTDBG_LEAK_CHECK_DF; // at end of program.
//      tmpFlag |= _CRTDBG_CHECK_ALWAYS_DF; // end of program
//      tmpFlag &= ~_CRTDBG_DELAY_FREE_MEM_DF; // turn it off

        //_CrtSetDbgFlag(_CRTDBG_ALLOC_MEM_DF | _CRTDBG_LEAK_CHECK_DF );

        //CRunDirector *rD = new COneGeneInhibitoryNetworkRunDirector();
        //CRunDirector *rD = new COneGeneInhibitoryAveragingDirector();
        //CRunDirector *rD = new CPC12EGFNGFRunDirector();
        //CEnsembleCombinationDirector *rD = new COneGeneEnsembleDirector(250);
        //CEnsembleCombinationDirector *rD = new CPC12EGFNGFEnsembleDirector(50);
        //CEnsembleCombinationDirector *rD = new CPC12AlternateEnsembleDirector(50);
//      CEnsembleCombinationDirector *rD = new CQuorumSensingEnsembleDirector(5);
//      CRunDirector *rD = new CQuorumSensingRunDirector();
        //CRunDirector *rD = new CLacRunDirector();
        //CRunDirector *rD = new CCombinatorialRunDirector();

//      rD->Execute();
//      delete rD;
//      cout << "Done plotting ensembledirector, proceed?" << endl;
//      cin.get();

    // define an uber-minimizable
        //NetworkMinimizableDirector *nM = new CPC12AlternateMinimizableDirector();

        //NetworkMinimizableDirector *nM = new PC12EGFNGFMinimizableDirector();

        //NetworkMinimizableDirector *nM = new CCoupledPC12NetworksMinimizableDirector(0.01,0.0);

//OneGeneInhibitoryNetwork is a quick test case for changes
        //NetworkMinimizableDirector *nM = new COneGeneInhibitoryNetworkMinimizableDirector(plotornot);
        //NetworkMinimizableDirector *nM = new COneGeneStochasticNetworkMinimizableDirector();

        NetworkMinimizableDirector *nM = new CQuorumSensingNetworkMinimizableDirector(plotornot);

        //NetworkMinimizableDirector *nM = new CCombinatorialNetworkMinimizableDirector();

        //NetworkMinimizableDirector *nM = new CAlternateOneGeneInhibitoryNetworkMinimizableDirector();

        //NetworkMinimizableDirector *nM = new CExampleNetworkMinimizableDirector();

        //NetworkMinimizableDirector *nM = new CCoupledOneGeneNetworksMinimizableDirector(0.01,0.0);

        //NetworkMinimizableDirector *nM = new CEGFRTraffickingCoupledMinimizableDirector(32.0,0.0);

        /*
        // read in two sets of parameters
        ParameterReader *pR1 = new ParameterReader("min2_cg_woPI3K.par");
        double *p1 = new double[nM->GetNParameters()];
        for(int i = 0; i < nM->GetNParameters(); i++)
        {
                p1[i] = pR1->ReadParameter();
        }
        delete pR1;
        ParameterReader *pR2 = new ParameterReader("min3_cg_woPI3K.par");
        double *p2 = new double[nM->GetNParameters()];
        for(i = 0; i < nM->GetNParameters(); i++)
        {
                p2[i] = pR2->ReadParameter();
        }
        delete pR2;
        ParameterReader *pR3 = new ParameterReader("min4_cg_woPI3K.par");
        double *p3 = new double[nM->GetNParameters()];
        for(i = 0; i < nM->GetNParameters(); i++)
        {
                p3[i] = pR3->ReadParameter();
        }
        delete pR3;
        
        // create the mapper and filter
        CLowDimensionalCostFunctionMapper *pMap = new CLowDimensionalCostFunctionMapper();
        CParameterFilter *pFilter = new CLogParameterFilter();
        // do the mapping
        pMap->Map(nM,pFilter,p1,p2,p3,-1,1,-1,1,961);

        delete pMap;
        delete pFilter;
        delete nM;
        delete [] p1;
        delete [] p2;
        delete [] p3;
        */

        double *parameters = new double[nM->GetNParameters()];
        cout << "Total number of parameters is " << nM->GetNParameters() << endl;
        if (parameterfile == defaultkey ) {
                // use default parameters
                for(i = 0; i < nM->GetNParameters(); i++)
                {
                        parameters[i] = nM->GetParameter(i);
                }
        }
        else {
                // read parameters from previous optimization pass
                ParameterReader *pReader = new ParameterReader(parameterfile);
                for(i = 0; i < nM->GetNParameters(); i++)
                {
                        parameters[i] = pReader->ReadParameter();
                }
                delete pReader;
        }

/*      // randomize parameters
        CParameterRandomizer *pR = new CParameterRandomizer(9891);
        pR->Randomize(parameters,nM->GetNParameters(),0.25,1);
*/
        //Print reactions in tex or mathematica input
/*      ReactionNetwork *rN = new CQuorumSensingNetwork();
        rN->GetReactions();
        rN->DumpTeXForm("QSequations.tex");
//      rN->DumpMathematicaForm("steadystate.dat");
        delete rN;
        cout << "texform printed\n";
        cin.get();
*/
        /*
        // randomize parameters
        CParameterRandomizer *pR = new CParameterRandomizer(776157);
        pR->Randomize(parameters,nM->GetNParameters(),0.5,3);
        */

        //cout << nM->ObjectiveFunction(parameters) << endl;

        /*
        // STEPSIZE TESTER - FOR USE WITH STOCHASTIC SENSITIVITY
        // evaluate the objective function at the minimum to determine what size step
        // to give to fixed stepsize mover
        double cost = nM->ObjectiveFunction(parameters);
        cout << "Cost of parameter set is " << cost << endl;
        for(int nEx = 0; nEx < nM->GetNExperiments(); nEx++)
        {
                cout << "Experiment " << nEx << ":" << endl;
                cout << "\tSmallest stepsize was " << nM->GetMover(nEx)->GetMinStepSize() << endl;
                cout << "\tLargest stepsize was " << nM->GetMover(nEx)->GetMaxStepSize() << endl;
                cout << "\tTotal number of steps was " << nM->GetMover(nEx)->GetTotalStepsTaken() << endl;
        }
        // just a dummy so I don't have to keep commenting out delete statement at bottom
        Minimizer *minimizer = new CBasicQuenchMinimizer(1,10,1.0,1);
        */
        
/*      CParameterFilter *pFilter = new CLogParameterFilter();
        CStochasticSensitivityAnalysis *minimizer = new CStochasticSensitivityAnalysis(false,true,25,2115,pFilter,1.0e-6);
        cout << "Performing Stochastic Sensitivity Analysis" << endl;
        minimizer->DoAnalysis(parameters,nM);
        cout << "RW Variance : " << minimizer->GetRWValue() << endl;
        delete pFilter;
*/
        /*
        cout << "(mu/rho,lambda)-ES One Selected" << endl;
        CParameterFilter *filter = new CLogParameterFilter();
        Minimizer *minimizer = new CCommaStrategyOne(filter,1.0,15,2,50,100,887);
        double cost = minimizer->Minimize(parameters,nM);
        cout << "Average cost of final population " << cost << endl;
        delete filter;
        */
        /*
        cout << "(mu,lambda)-ES Selected" << endl;
        CParameterFilter *filter = new CLogParameterFilter();
        Minimizer *minimizer = new CClonalCommaES(filter,0.0,15,100,1109);
        double cost = minimizer->Minimize(parameters,nM);
        cout << "Average cost of final population " << cost << endl;
        delete filter;
        */

        // BE CAREFUL WITH THE ERROR TOLERANCE, PARTICULARLY IF YOU USE A LOW ORDER
        // INTEGRATOR
//      CParameterFilter *filter = new CIdentityParameterFilter();
        CParameterFilter *filter = new CLogParameterFilter();
        double cost;
        Minimizer *minimizer;
        int nEquilibrationSteps;
        double stepScale;
        CLeastSquaresObserver *LMForcesObserver = new CLeastSquaresObserver();
        switch (smethod) {
        case 'l':
                cout << "Levenberg-Marquardt Selected" << endl;
                minimizer = new CImprovedLevenbergMarquardtMinimizer(filter,false,1.0e-06,false,1.0e-08,true,1.0e-06,true,1.0e-06,1000,1.0e-06);
                //CLeastSquaresObserver *LMForcesObserver = new CLeastSquaresObserver();
                minimizer->Attach(LMForcesObserver);
                cost = minimizer->Minimize(parameters,nM);
                cout << "Final cost of " << cost << endl;
                minimizer->Detach(LMForcesObserver);
                delete filter;
                break;

        case 's':
                cout << "Scale-Invariant Levenberg-Marquardt Selected" << endl;
                minimizer = new CScaleInvariantMarquardtMinimizer(filter,false,1.0e-2,true,1.0e-08,true,1.0e-06,true,1.0e-06,1000,1.0e-06);
                cost = minimizer->Minimize(parameters,nM);
                cout << "Final cost of " << cost << endl;
                delete filter;
                break;

        case 'c':
                cout << "Conjugate Gradient Selected" << endl;
                minimizer = new CConjugateGradientMinimizer(false,filter,3,1.0e-06);
//      CConjugateGradientObserver *CGForcesObserver = new CConjugateGradientObserver();
//      minimizer->Attach(CGForcesObserver);
                cost = minimizer->Minimize(parameters,nM);
//      minimizer->Detach(CGForcesObserver);
                delete filter;
                break;

        case 'p':
//      GnuPlotterSimulatedAnnealingObserver *costObserver = new GnuPlotterSimulatedAnnealingObserver();
//      costObserver->SetTitle("Basic Quench");
                cout << "Plot Selected" << endl;
                nEquilibrationSteps = 25;
                stepScale = 0.0;
                minimizer = new CBasicQuenchMinimizer(nEquilibrationSteps,10,stepScale,99449,filter);
//      minimizer->Attach(costObserver);
                cost = minimizer->Minimize(parameters,nM);
//      minimizer->Detach(costObserver);
//      delete costObserver;
                break;

        case 'b':
//      GnuPlotterSimulatedAnnealingObserver *costObserver = new GnuPlotterSimulatedAnnealingObserver();
//      costObserver->SetTitle("Basic Quench");
                cout << "Basic Quench Selected" << endl;
                nEquilibrationSteps = 25;
                stepScale = 1.0;
                minimizer = new CBasicQuenchMinimizer(nEquilibrationSteps,10,stepScale,99449,filter);
//      minimizer->Attach(costObserver);
                cost = minimizer->Minimize(parameters,nM);
//      minimizer->Detach(costObserver);
//      delete costObserver;
                break;

        case 'h':
                CParameterFilter *pFilter = new CLogParameterFilter();
                CStochasticSensitivityAnalysis *minimizer = new CStochasticSensitivityAnalysis(false,true,50,2115,pFilter,1.0e-6);
                cout << "Performing Stochastic Sensitivity Analysis" << endl;
                minimizer->DoAnalysis(parameters,nM);
                cout << "RW Variance : " << minimizer->GetRWValue() << endl;
                delete pFilter;
                break;

        }//end switch method

        cout << "Best cost for this minimization is " << cost << endl;
        // write the best parameters to a file and stdout
        fstream *pfbest = new fstream("goodtemp10.par",ios::out);
        *pfbest << "##################################################################" << endl;
        *pfbest << "# BEST PARAMETERS" << endl;
        *pfbest << "# COST : " << cost << endl;
        *pfbest << "##################################################################" << endl;
        for(i = 0; i < nM->GetNParameters(); i++)
        {
                *pfbest << parameters[i] << endl;
        }
        delete pfbest;

        cout << "finish?" << endl;
        cin.get();

        delete minimizer;
        delete [] parameters;
        delete nM;

        return 0;

}

---