#!/bin/bash

module load xfdtd #load the XF module for use

cd Evolved_Dipole #move to Evolved_Dipole Directory

./Evolved_Dipole --start #Create First Gen

while inotifywait ~/Evolved_Dipole/watches -r -e modify -m; do #wait until watches directory is modified
    if tail -n1 ~/Evolved_Dipole/watches/watch.txt | grep 0; then #if a 0 is written to watch.txt then run xf macro
	xfui --execute-macro-script=/users/PAS0654/osu8742/XFScripts/dipole_PEC.xmacro
    elif tail -n1 ~/Evolved_Dipole/watches/watch.txt | grep 1; then #if a 1 has been written to watch.txt then make Gen 2 or higher
	./Evolved_Dipole --cont
    else #if anything else is written then end
    fi
done

//Macro to build a dipole antenna made of PEC in XF
//J. Potter -- 06/17

var path = "/users/PAS0654/osu8742/Evolved_Dipole/handshake.csv";
var file = new File(path);
file.open(1);
var hands = file.readAll();
//hands=hands.replace(/,/g, '');
//hands=hands.replace(/[(]/g, '');
//hands=hands.replace(/[)]/g, '');
//Output.println(hands);
var radii=[];
var lengths=[];


//Only take the lines with functions
var lines = hands.split('\n');
for(var i = 0;i < lines.length;i++){
    if(i>=8)
    {
	var params = lines[i].split(",");
	    radii[i-8]=params[0];
		lengths[i-8]=params[1];
			Output.println(radii[i-8]);
				Output.println(lengths[i-8]);
					
					}
}
file.close();

// set variables for length of the dipoles, radius of the dipoles and units for the lengths
var dipoleLength = 10;
var dipoleRadius = 5;
var connectorLength=1.0;
var gridSize = .25;
var units = " cm" //don't delete the space
var frequency = .5;

//App.saveCurrentProjectAs("C:/Users/Jordan/Desktop/Emacs_Scripts/XFScripts/XFDipoleTest1.xf");


function CreateDipole()
{
	App.getActiveProject().getGeometryAssembly().clear();

	//create a new sketch
	var dipole = new Sketch();
	var base = new Ellipse( new Cartesian3D(0,0,0), new Cartesian3D( dipoleRadius + units,0,0 ), 1.0, 0.0, Math.PI*2 );
	dipole.addEdge(base);

	//add depth to the circle
	var extrude = new Extrude( dipole, dipoleLength + units );

	//create a recipe and model -- still need to figure what a recipe is... 
	var dipoleRecipe = new Recipe();
	dipoleRecipe.append(extrude);
	var dipoleModel = new Model();
	dipoleModel.setRecipe(dipoleRecipe);
	dipoleModel.name = "Dipole Antenna Test";

	//set locations of the left and right dipole segments
	dipoleModel.getCoordinateSystem().translate( new Cartesian3D(0,0,1 + units));
	var dipoleInProject1 = App.getActiveProject().getGeometryAssembly().append(dipoleModel);
	dipoleModel.getCoordinateSystem().translate(new Cartesian3D(0,0,(-1 - dipoleLength) + units));
	var dipoleInProject2 = App.getActiveProject().getGeometryAssembly().append(dipoleModel);

	// Now set the material for the Dipole:
    	var pecMaterial = App.getActiveProject().getMaterialList().getMaterial( "PEC" );
    	if( null == pecMaterial )
    	{
		Output.println( "\"PEC\" material was not found, could not associate with the antenna." );
    	}
    	else
    	{
		App.getActiveProject().setMaterial( dipoleInProject1, pecMaterial );
		App.getActiveProject().setMaterial( dipoleInProject2, pecMaterial );
    	}

	//zoom to view the extent of the creation
	View.zoomToExtents();
}

function CreatePEC() //borrowed from XF demo
{
    //Make the material.  We will use PEC, or Perfect Electrical Conductor:
    var pec = new Material();
    pec.name = "PEC";
    
    var pecProperties = new PEC();      // This is the electric properties that defines PEC
    var pecMagneticFreespace = new MagneticFreespace();     // We could make a material that acts as PEC and PMC, but in this case we just care about electrical components.
    var pecPhysicalMaterial = new PhysicalMaterial();
    pecPhysicalMaterial.setElectricProperties( pecProperties );
    pecPhysicalMaterial.setMagneticProperties( pecMagneticFreespace );
    pec.setDetails( pecPhysicalMaterial );

    // PEC is historically a "white" material, so we can easily change its appearance:
    var pecBodyAppearance = pec.getAppearance();
    var pecFaceAppearance = pecBodyAppearance.getFaceAppearance();  // The "face" appearance is the color/style associated with the surface of geometry objects
    pecFaceAppearance.setColor( new Color( 255, 255, 255, 255 ) );  // Set the surface color to white. (255 is the maximum intensity, these are in order R,G,B,A).
    
    // Check for an existing material
    if( null != App.getActiveProject().getMaterialList().getMaterial( pec.name ) )
    {
        App.getActiveProject().getMaterialList().removeMaterial( pec.name );
    }
    App.getActiveProject().getMaterialList().addMaterial( pec );
}


function CreateAntennaSource()
{
    // Here we will create our waveform, create our circuit component definition for the feed, and create
    // a CircuitComponent that will attach those to our current geometry.
    var waveformList = App.getActiveProject().getWaveformList();
    // clear the waveform list
    waveformList.clear();

    var parameterList = App.getActiveProject().getParameterList();
    parameterList.clear();
    parameterList.addParameter("freq",".5");
    
    // Create a sinusoidal input wave
    var waveform = new Waveform();
    var sine = new RampedSinusoidWaveformShape();
    //sine.setFrequency("freq" + " GHz");  //Uncomment if using parameter sweep
    sine.setFrequency(frequency + " GHz"); //Comment if no parameter sweep
    waveform.setWaveformShape( sine );
    waveform.name ="Sinusoid";
    var waveformInList = waveformList.addWaveform( waveform );

    // Now to create the circuit component definition:
    var componentDefinitionList = App.getActiveProject().getCircuitComponentDefinitionList();
    // clear the list
    componentDefinitionList.clear();

    // Create our Feed
    var feed = new Feed();
    feed.feedType = Feed.Voltage; // Set its type enumeration to be Voltage.
    // Define a 50-ohm resistance for this feed
    var rlc = new RLCSpecification();
    rlc.setResistance( "50 ohm" );
    rlc.setCapacitance( "0" );
    rlc.setInductance( "0" );
    feed.setImpedanceSpecification( rlc );
    feed.setWaveform( waveformInList );  // Make sure to use the reference that was returned by the list, or query the list directly
    feed.name = "50-Ohm Voltage Source";
    var feedInList = componentDefinitionList.addCircuitComponentDefinition( feed );

    // Now create a circuit component that will be the feed point for our simulation
    var componentList = App.getActiveProject().getCircuitComponentList();
    componentList.clear();
    
    var component = new CircuitComponent();
    component.name = "Source";
    component.setAsPort( true );
    // Define the endpoints of this feed - these are defined in world position, but you can also attach them to edges, faces, etc.
    var coordinate1 = new CoordinateSystemPosition( 0, 0,0 );
    var coordinate2 = new CoordinateSystemPosition( 0, 0, "1" + units );
    component.setCircuitComponentDefinition( feedInList );
    component.setEndpoint1( coordinate1 );
    component.setEndpoint2( coordinate2 );
    componentList.addCircuitComponent( component );
}

function CreateGrid()
{
    // Set up the grid spacing for the dipole antenna
    var grid = App.getActiveProject().getGrid();
    var cellSizes = grid.getCellSizesSpecification();
    

    cellSizes.setTargetSizes( Cartesian3D( gridSize + units, gridSize + units, gridSize + units ) );
    // And we need to set the Minimum Sizes - these are the minimum deltas that we will allow in this project.
    // We'll use the scalar ratio of 20% here.
    cellSizes.setMinimumSizes( Cartesian3D( ".5", ".5", ".5" ) );
    cellSizes.setMinimumIsRatioX( true );
    cellSizes.setMinimumIsRatioY( true );
    cellSizes.setMinimumIsRatioZ( true );

    grid.specifyPaddingExtent( Cartesian3D( "20", "20", "20" ), Cartesian3D( "20", "20", "20" ), true, true );
}

function CreateSensors()
{
    // Here we will create a sensor definition and attach it to a near-field sensor on the surface of one of our objects.

    var sensorDataDefinitionList = App.getActiveProject().getSensorDataDefinitionList();
    sensorDataDefinitionList.clear();

    // Create a sensor
    var farSensor = new FarZoneSensor();
    farSensor.retrieveTransientData = true;
    farSensor.setAngle1IncrementRadians(Math.PI/12.0);
    farSensor.setAngle2IncrementRadians(Math.PI/12.0);
    farSensor.name = "Far Zone Sensor";
 

    var FarZoneSensorList = App.getActiveProject().getFarZoneSensorList();
    FarZoneSensorList.clear();
    FarZoneSensorList.addFarZoneSensor( farSensor );
}



function CreateAntennaSimulationData()
{
    // This function modifies the NewSimulationData parameters of the project.
    // They're called "New" because they get applied every time we create an instance of this
    // project and place it on our simulation queue.
    var simData = App.getActiveProject().getNewSimulationData();
    
    // These should already be set, however just to make sure the current project is set up correctly
    simData.excitationType = NewSimulationData.DiscreteSources;
    simData.enableSParameters = false;

    // Set convergence threshold to -40 dB, and maximum number of timesteps to 10000.
    var terminationCriteria = simData.getTerminationCriteria();
    terminationCriteria.convergenceThreshold = -40;
    terminationCriteria.setMaximumSimulationTime( "20000*timestep" );

    // Do not attempt to collect steady-state data
    simData.getFOIParameters().collectSteadyStateData = true;
 
    // Construct parameter sweep
    //var sweep = simData.getParameterSweep();
    //sweep.parameterName = "freq";
    //sweep.setParameterValuesByCount(.1, 2, 5); //(start value ,end value, # of steps)
    //simData.enableParameterSweep = true;


}

function QueueSimulation()
{
    // Create and start the simulation. Project must be saved in order to run.
    var simulation = App.getActiveProject().createSimulation( true );   

    Output.println( "Successfully created the simualtion." );

    var projectId = simulation.getProjectId();
    var simulationId = simulation.getSimulationId();
    var numRuns = simulation.getRunCount();
}

for(var i = 0;i < 5;i++){
	var dipoleLength = lengths[i];
	var dipoleRadius = radii[i];
//Actually create the material and then the dipole
	CreatePEC();
	CreateDipole();
	CreateAntennaSource();
	CreateGrid();
	CreateSensors();
	CreateAntennaSimulationData();
	QueueSimulation();
}

/* 
 * File:   Evolved_Dipole.cpp
 * Author: Jordan Potter
 *
 * Created on July 14, 2017, 11:07 AM
 */

#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include <string>
#include <sstream>

using namespace std;

const int NVAR=2; //number of variables (r and l)
const int NPOP=5; //size of population


double rand_normal(double mean, double stddev)
{//Box muller method to generate normal numbers
    static double n2 = 0.0;
    static int n2_cached = 0;
    if (!n2_cached)
    {
        double x, y, r;
        do
        {
            x = 2.0*rand()/RAND_MAX - 1;
            y = 2.0*rand()/RAND_MAX - 1;

            r = x*x + y*y;
        }
        while (r == 0.0 || r > 1.0);
        {
            double d = sqrt(-2.0*log(r)/r);
            double n1 = x*d;
            n2 = y*d;
            double result = n1*stddev + mean;
            n2_cached = 1;
            return result;
        }
    }
    else
    {
        n2_cached = 0;
        return n2*stddev + mean;
    }
}



void GenerateSample(double new_val[][NVAR], double mean_val[], double std_dev[],int n){
    double u;
    for(int i=0;i<NVAR;i++)
    {
        //Create a sample of mean + stdev*(gaussian sample) to use
        u=rand_normal(0,.5);
        new_val[n][i]=mean_val[i]+std_dev[i]*u;
        
        //Include conditionals for bounds and constraints here
        int errors=0;
        if(new_val[n][i] <= 0)
        {
            i--; //If a parameter is less than 0, decrement so this value is created again
            errors++;
        }
        
        if(errors >= 50) //If there are 50 errors something is going wrong, so stop
        {
            cout << "Error: Sample Generation Failed" << endl;
            return;
        }
    }
}
void Simulation(double x[][NVAR], double y[NVAR], double mean[], double deviation[], double best[]){
    //Write to CSV file. Might need to slightly adjust for generations after #1...
    //For now, this will only write the last generation to file
    ofstream handshake;
    handshake.open("handshake.csv");
    handshake << "C++ ESTRA -- Jordan Potter" << "\n";
    handshake << "Mean Vector:" << "\n"; //Write the Current Mean Values to handshake.csv
    for(int i=0;i<NVAR;i++)
      {
	if(i==(NVAR-1))
	  {
	    handshake << mean[i] << "\n";
	  }
	else
	  {
	    handshake << mean[i] << ",";
	  }
      }
    handshake << "Deviation Vector:" << "\n"; //Write the Current Deviation Values to handshake.csv
    for(int i=0;i<NVAR;i++)
      {
	if(i==(NVAR-1))
	  {
	    handshake << deviation[i] << "\n";
	  }
	else
	  {
	    handshake << deviation[i] << ",";
	  }
      }
    handshake << "Best Antenna Vector:" << "\n"; //Write the Current Deviation Values to handshake.csv
    for(int i=0;i<(NVAR+1);i++)
      {
	if(i==(NVAR))
	  {
	    handshake << best[i] << "\n";
	  }
	else
	  {
	    handshake << best[i] << ",";
	  }
      }
    handshake << "Generation:" << "\n"; //Write the Generation Values to handshake.csv
    for(int i=0;i<NPOP;i++)
    {
        for(int j=0;j<NVAR;j++)
        {
            if(j==(NVAR-1))
            {
                handshake << x[i][j] << "\n";
            }
            else
            {
                handshake << x[i][j] << ",";
            }
        }
    }
    handshake.close();

    //Set Value in another file, so that controller.sh knows to start XF
}

void Read(double x[][NVAR], double y[NVAR], double mean[], double deviations[], double best[]){    
    //Import values from Handshake for mean/x population/current deviations
  ifstream handshake;
  handshake.open("handshake.csv");
  string csvContent[NPOP+8]; //contain each line of csv
  string strToDbl; //gets overwritten to contain each cell of a line. Then transferred to x,y,mean,deviations,best
  for(int i=0;i<(NPOP+8);i++)
    {
      getline(handshake,csvContent[i]); //read in mean
      if(i==2)
	{
	  istringstream stream(csvContent[i]);
	  for(int j=0;j<NVAR;j++)
	    {
	      getline(stream,strToDbl,',');
	      mean[j] = atof(strToDbl.c_str());
	    }
	}
      else if(i==4) //read in deviations
	{
	  istringstream stream(csvContent[i]);
	  for(int j=0;j<NVAR;j++)
	    {
	      getline(stream,strToDbl,',');
	      deviations[j] = atof(strToDbl.c_str());
	    }
	}
      else if(i==6) //read in the current best values
	{
	  istringstream stream(csvContent[i]);
	  for(int j=0;j<NVAR+1;j++)
	    {
	      getline(stream,strToDbl,',');
	      best[j] = atof(strToDbl.c_str());
	    }
	}
      else if(i>=8) //read in the generation
	{
	  istringstream stream(csvContent[i]);
	  for(int j=0;j<NVAR;j++)
	    {
	      getline(stream,strToDbl,',');
	      x[i-8][j] = atof(strToDbl.c_str());
	    }
	}
    }


    //Import Values from Simulation from Handshook that XF (or eventually ARASim) will write
    //May need to adjust if format of handshook changes
  ifstream handshook;
  handshook.open("handshook.csv");
  string strToDbl2[NPOP+2];
  for(int i=0;i<(NPOP+2);i++)
    {
      getline(handshook,strToDbl2[i]);
      if(i>=2)
        {
	  y[i-2] = atof(strToDbl2[i].c_str());
        }
    }
  
  handshook.close();
}

int Mutate(double x[][NVAR],double out[], double best[], double mean[]){
    //check to see if the output from x[i][] is larger than the output from x[m].
    //i.e. output[i]>mean_output
    int count=0;
    for(int i=0;i<NPOP;i++)
    {
        if(out[i]>best[0])
        {
            best[0]=out[i];
            for(int j=0;j<NVAR;j++)
            {
                best[j+1]=x[i][j];
                mean[j]=x[i][j];
            }
            count++;
        }
    }
    return count;
}
void CheckConvergence(int numSuccess, double p, double q, double d[]){
    //If things aren't getting better than expand search radius. 
    //However, if things are getting better tighten search radius.
    //Check this over k? values of x.
    //P(output)>p then d=d/q
    //if not then d=d*q
    double P = (double) numSuccess/NPOP;
    if(P >= p)
    {
        for(int i=0;i<NVAR;i++)
        {
            d[i]=d[i]/q;
        }
    }
    else
    {
        for(int i=0;i<NVAR;i++)
        {
            d[i]=d[i]*q;
        }
    }
}


int main(int argc, char** argv) {
    double m[NVAR]={.5,.5}; //mean values of r and l
    double x[NPOP][NVAR]; //produced value of r and l
    double d[NVAR] = {.25,.25}; //standard deviation vector
    double best[NVAR+1]={-40,0,0}; //array to store param values and output score of top antenna {fit_score,r,l}
    double output[NPOP]; //array to store scores of each generations
    double q = .9; //factor to adjust search radius by
    double p = .2; //ratio of samples that must be correct to adjust search radius
    int numSuccess; //number of samples that improve on the best value
    srand(time(0));
    
    if(argc != 2)
        cout << "Error: Usage. Specify --start or --cont" << endl;
    else
    {
      if(string(argv[1]) == "--start")
      {

	for(int i=0;i<NPOP;i++)
	  {
	    GenerateSample(x,m,d,i); //Generate NPOP samples
	  }
	Simulation(x, output, m, d, best); // Write current population to disk as well as mean,deviations and best population
      }
      else if(string(argv[1]) == "--cont")
        {
	  //Read Function (read in population to x[],mean to m[],deviations to d[], best to best[] and scores to output[]. 
	  Read(x,output,m,d,best);
	  numSuccess = Mutate(x, output, best, m);
	  CheckConvergence(numSuccess, p, q, d);
	  

	  //Check if d/d_o<lambda has been achieved. Consider adding this to check convergence function. 
	  
	  cout << "The best parameter values are:" << endl;
	  for(int i =0; i <NVAR; i++)
	    {
	      cout << m[i] << endl; //Consider writing to a log file
            }

	  for(int i=0;i<NPOP;i++)
	    {
	      GenerateSample(x,m,d,i);
	    }
	  Simulation(x, output, m, d, best);
        }
      else
        {
	  cout << "Error: Specify --start or --cont" << endl;
        }
    }
    return 0;
}

var query = new ResultQuery(); 
///////////////////////Get Theta and Phi Gain///////////////
query.projectId = App.getActiveProject().getProjectDirectory(); 
query.simulationId = "000062";
query.runId = "Run0001"; 
query.sensorType = ResultQuery.FarZoneSensor;
query.sensorId = "Far Zone Sensor"; 
query.timeDependence = ResultQuery.SteadyState; 
query.resultType = ResultQuery.Gain; 
query.fieldScatter = ResultQuery.TotalField;
query.resultComponent = ResultQuery.Theta;
query.dataTransform = ResultQuery.NoTransform;
query.complexPart = ResultQuery.NotComplex;
query.surfaceInterpolationResolution = ResultQuery.NoInterpolation;
query.setDimensionRange( "Frequency", 0, -1 );
query.setDimensionRange( "Theta", 0, -1 );
query.setDimensionRange( "Phi", 0, -1 );

var thdata = new ResultDataSet( "" ); 
thdata.setQuery( query ); 
if( !thdata.isValid() ){ 
Output.println( "1getCurrentDataSet() : " + 
thdata.getReasonWhyInvalid() ); 
}

query.resultComponent = ResultQuery.Phi;
var phdata = new ResultDataSet("");
phdata.setQuery(query);

if( !phdata.isValid() ){ 
Output.println( "2getCurrentDataSet() : " + 
phdata.getReasonWhyInvalid() ); 
}
/////////////////Get Theta and Phi Phase///////////////////////////////////

query.resultType = ResultQuery.E; 
query.fieldScatter = ResultQuery.TotalField;
query.resultComponent = ResultQuery.Theta;
query.dataTransform = ResultQuery.NoTransform;
query.complexPart = ResultQuery.Phase;
query.surfaceInterpolationResolution = ResultQuery.NoInterpolation;
query.setDimensionRange( "Frequency", 0, -1 );
query.setDimensionRange( "Theta", 0, -1 );
query.setDimensionRange( "Phi", 0, -1 );


var thphase = new ResultDataSet("");
thphase.setQuery(query);

if( !thphase.isValid() ){ 
Output.println( "3getCurrentDataSet() : " + 
thphase.getReasonWhyInvalid() ); 
}

query.resultComponent = ResultQuery.Phi;
query.ComplexPart = ResultQuery.Phase;
var phphase = new ResultDataSet("");
phphase.setQuery(query);

if( !phphase.isValid() ){ 
Output.println( "4getCurrentDataSet() : " + 
phphase.getReasonWhyInvalid() ); 
}

/////////////////Get Input Power///////////////////////////
query.sensorType = ResultQuery.System;
query.sensorId = "System"; 
query.timeDependence = ResultQuery.SteadyState; 
query.resultType = ResultQuery.NetInputPower; 
query.fieldScatter = ResultQuery.NoFieldScatter;
query.resultComponent = ResultQuery.Scalar;
query.dataTransform = ResultQuery.NoTransform;
query.complexPart = ResultQuery.NotComplex;
query.surfaceInterpolationResolution = ResultQuery.NoInterpolation;
query.clearDimensions();
query.setDimensionRange("Frequency",0,-1);

var inputpower = new ResultDataSet("");
inputpower.setQuery(query);
if( !inputpower.isValid() ){ 
Output.println( "5getCurrentDataSet() : " + 
inputpower.getReasonWhyInvalid() ); 
}

FarZoneUtils.exportToUANFile(thdata,thphase,phdata,phphase,inputpower,"/users/PAS0654/osu8742/Evolved_Dipole/data/1.uan");