"

From 2010.igem.org

1. include <iostream>
2. include <vector>
3. include <fstream>
4. include <math.h>

1. ifndef __mjdmatrix_h
2. define __mjdmatrix_h
3. include <iostream>

using namespace std; // generic object (class) definition of matrix: template <class D> class matrix{

 // NOTE: maxsize determines available memory storage, but
 // actualsize determines the actual size of the stored matrix in use
 // at a particular time.
 int maxsize;  // max number of rows (same as max number of columns)
 int actualsize;  // actual size (rows, or columns) of the stored matrix
 D* data;      // where the data contents of the matrix are stored
 void allocate()   {
   delete[] data;
   data = new D [maxsize*maxsize];
   };
 matrix() {};                  // private ctor's
 matrix(int newmaxsize) {matrix(newmaxsize,newmaxsize);};

public:

void print() { for(int i=0; i<maxsize; i++){ for(int j=0; j<maxsize; j++){

cout << getvalue(i,j) << " "; } cout << endl; }

};

D getvalue(int row, int column){return data[ row * maxsize + column ];};

 matrix(int newmaxsize, int newactualsize)  { // the only public ctor
   if (newmaxsize <= 0) newmaxsize = 5;
   maxsize = newmaxsize; 
   if ((newactualsize <= newmaxsize)&&(newactualsize>0))
     actualsize = newactualsize;
   else 
     actualsize = newmaxsize;
   // since allocate() will first call delete[] on data:
   data = 0;
   allocate();
   };
 ~matrix() { delete[] data; };

 void comparetoidentity()  {
   int worstdiagonal = 0;
   D maxunitydeviation = 0.0;
   D currentunitydeviation;
   for ( int i = 0; i < actualsize; i++ )  {
     currentunitydeviation = data[i*maxsize+i] - 1.;
     if ( currentunitydeviation < 0.0) currentunitydeviation *= -1.;
     if ( currentunitydeviation > maxunitydeviation )  {
          maxunitydeviation = currentunitydeviation;
          worstdiagonal = i;
          }
     }
   int worstoffdiagonalrow = 0;
   int worstoffdiagonalcolumn = 0;
   D maxzerodeviation = 0.0;
   D currentzerodeviation ;
   for ( int i = 0; i < actualsize; i++ )  {
     for ( int j = 0; j < actualsize; j++ )  {
       if ( i == j ) continue;  // we look only at non-diagonal terms
       currentzerodeviation = data[i*maxsize+j];
       if ( currentzerodeviation < 0.0) currentzerodeviation *= -1.0;
       if ( currentzerodeviation > maxzerodeviation )  {
         maxzerodeviation = currentzerodeviation;
         worstoffdiagonalrow = i;
         worstoffdiagonalcolumn = j;
         }

       }
     }
   cout << "Worst diagonal value deviation from unity: " 
      << maxunitydeviation << " at row/column " << worstdiagonal << endl;
   cout << "Worst off-diagonal value deviation from zero: " 
      << maxzerodeviation << " at row = " << worstoffdiagonalrow 
      << ", column = " << worstoffdiagonalcolumn << endl;
 }

 void settoproduct(matrix& left, matrix& right)  {
   actualsize = left.getactualsize();
   if ( maxsize < left.getactualsize() )   {
     maxsize = left.getactualsize();
     allocate();
     }
   for ( int i = 0; i < actualsize; i++ )
     for ( int j = 0; j < actualsize; j++ )  {
       D sum = 0.0;
       D leftvalue, rightvalue;
       bool success;
       for (int c = 0; c < actualsize; c++)  {
           left.getvalue(i,c,leftvalue,success);
           right.getvalue(c,j,rightvalue,success);
           sum += leftvalue * rightvalue;
           }
       setvalue(i,j,sum);
       }
   }

 void copymatrix(matrix&  source)  {
   actualsize = source.getactualsize();
   if ( maxsize < source.getactualsize() )  {
     maxsize = source.getactualsize();
     allocate();
     }
   for ( int i = 0; i < actualsize; i++ )
     for ( int j = 0; j < actualsize; j++ )  {
       D value;
       bool success;
       source.getvalue(i,j,value,success);
       data[i*maxsize+j] = value;
       }
   };

 void setactualsize(int newactualsize) {
   if ( newactualsize > maxsize )
     {
       maxsize = newactualsize ; // * 2;  // wastes memory but saves
                                     // time otherwise required for
                                     // operation new[]
       allocate();
     }
   if (newactualsize >= 0) actualsize = newactualsize;
   };

 int getactualsize() { return actualsize; };

 void getvalue(int row, int column, D& returnvalue, bool& success)   {
   if ( (row>=maxsize) || (column>=maxsize) 
     || (row<0) || (column<0) )
     {  success = false;
        return;    }
   returnvalue = data[ row * maxsize + column ];
   success = true;
   };

 bool setvalue(int row, int column, D newvalue)  {
   if ( (row >= maxsize) || (column >= maxsize) 
       || (row<0) || (column<0) ) return false;
   data[ row * maxsize + column ] = newvalue;
   return true;
   };

 void invert()  {
   if (actualsize <= 0) return;  // sanity check
   if (actualsize == 1) return;  // must be of dimension >= 2
   for (int i=1; i < actualsize; i++) data[i] /= data[0]; // normalize row 0
   for (int i=1; i < actualsize; i++)  { 
     for (int j=i; j < actualsize; j++)  { // do a column of L
       D sum = 0.0;
       for (int k = 0; k < i; k++)  
           sum += data[j*maxsize+k] * data[k*maxsize+i];
       data[j*maxsize+i] -= sum;
       }
     if (i == actualsize-1) continue;
     for (int j=i+1; j < actualsize; j++)  {  // do a row of U
       D sum = 0.0;
       for (int k = 0; k < i; k++)
           sum += data[i*maxsize+k]*data[k*maxsize+j];
       data[i*maxsize+j] = 
          (data[i*maxsize+j]-sum) / data[i*maxsize+i];
       }
     }
   for ( int i = 0; i < actualsize; i++ )  // invert L
     for ( int j = i; j < actualsize; j++ )  {
       D x = 1.0;
       if ( i != j ) {
         x = 0.0;
         for ( int k = i; k < j; k++ ) 
             x -= data[j*maxsize+k]*data[k*maxsize+i];
         }
       data[j*maxsize+i] = x / data[j*maxsize+j];
       }
   for ( int i = 0; i < actualsize; i++ )   // invert U
     for ( int j = i; j < actualsize; j++ )  {
       if ( i == j ) continue;
       D sum = 0.0;
       for ( int k = i; k < j; k++ )
           sum += data[k*maxsize+j]*( (i==k) ? 1.0 : data[i*maxsize+k] );
       data[i*maxsize+j] = -sum;
       }
   for ( int i = 0; i < actualsize; i++ )   // final inversion
     for ( int j = 0; j < actualsize; j++ )  {
       D sum = 0.0;
       for ( int k = ((i>j)?i:j); k < actualsize; k++ )  
           sum += ((j==k)?1.0:data[j*maxsize+k])*data[k*maxsize+i];
       data[j*maxsize+i] = sum;
       }
   };

};

1. endif

using namespace std;

class Function { private: vector <double> wavelength; vector <double> function; vector <double> lambda;

double RBF(double, double); double sigma; double width; int SIZE;

public: void getData(const char*); void findLambda(double); void fix();

double size(){return wavelength.size();} double lat(int i){return lambda.at(i);} double wat(int i){return wavelength.at(i);} double fat(int i){return function.at(i);} void wequals(int i, double value){wavelength.at(i)=value;} void fequals(int i, double value){function.at(i)=value;}

void eraseBeginning(int); void erase(int);

void print(); void print(double,double,double); void printStore();

double value(double);

};

double X(double j, double a, double h) {return a+j*h;} double H(double a,double b,double n) {return (b-a)/n;}

double integral (double a, double b, Function);

void runProgram(char*,char*,double,double);

double luminousFlux(Function, Function,double,double);

void sortData(char*, char*); void testData(char* data,double, char* output);

void brightness(double lux);

1. define MINI 200
2. define MAXI 800
3. define HCON 6.63e-34
4. define CCON 2.99e8
5. define TURNOVER 0.125*2000
6. define CDEN 8e12

1. ifndef FILE

#define FILE "foo.dat"

1. endif

1. define N 1000

int main() { //sortData("workingdata.dat","spectra.dat"); //testData("photopic.dat",0.7,"output3.dat"); runProgram("scotopic.dat","spectra.dat",CDEN,TURNOVER); return 0; }

void brightness(double lux) { if(lux<10e-5) cout << "invisible" << endl; else if(lux<10e-4) cout << "the light from Sirius, the brightest star in the night sky" << endl; else if(lux<0.002) cout << "total starlight with an overcast sky" << endl; else if(lux<0.01) cout << "a moonless clear night sky with airglow" << endl; else if(lux<0.27) cout << "a night sky with a quarter moon" << endl; else if(lux<1) cout << "a full moon on a clear night" << endl; else if(lux<3.4) cout << "a full moon overhead at tropical latitude" << endl; else if(lux<50) cout << "the dark limit of civil twilight under a clear sky" << endl; else if(lux<80) cout << "a family living room" << endl; else if(lux<100) cout << "a hallway/toilet" << endl; else if(lux<320) cout << "a very dark overcast day" << endl; else if(lux<1e3) cout << "office lighting or a sunrise/sunset" << endl; else if(lux<1e4) cout << "an overcast day" << endl; else if(lux<32e3) cout << "full daylight (not direct sun)" << endl; else cout << "Direct sunlight " << endl;

}

void sortData(char* input, char* output) { Function a; a.getData(input); double wMIN=a.wat(0), wMAX=a.fat(0), fMIN=a.wat(1), fMAX=a.fat(1); a.eraseBeginning(4); double xorigin=(a.wat(0)+a.wat(2))/2, yorigin=(a.fat(0)+a.fat(1))/2; double xMAX=a.wat(1), yMAX=a.fat(2); a.eraseBeginning(6);

double yratio= (fMAX-fMIN)/(yMAX-yorigin); double xratio= (wMAX-wMIN)/(xMAX-xorigin);

for(int i=0; i<a.size(); i++){ a.wequals(i, wMIN+(a.wat(i)-xorigin)*xratio); a.fequals(i, fMIN+(a.fat(i)-yorigin)*yratio); }

int i; a.fix();

ofstream fout(output);

if(!fout) {cout << "Could not open file " << output << " program terminated." << endl; fout.close(); return;}

for(double i=0; i<a.size() ; i++) fout << a.wat(i) << endl << a.fat(i) << endl;

fout.close();

}

void Function::eraseBeginning(int i) { i/=2; function.erase(function.begin(), function.begin()+i); wavelength.erase(wavelength.begin(),wavelength.begin()+i); }

void runProgram(char* lumFile, char* specFile,double cden, double turnover ) { Function luminosity, spectra;

cout << "Getting luminosity data..." << endl; luminosity.getData(lumFile); cout << "Calculating luminosity function..." << endl; luminosity.findLambda(0.7);

cout << "Getting spectral data..." << endl; spectra.getData(specFile); cout << "Calculating spectral function..." << endl; spectra.findLambda(0.7);

cout << "Calculating luminous flux..." << endl; double x=luminousFlux(luminosity,spectra,cden,turnover) ;

cout << "Luminous Flux is: " << x << endl << "This would be as bright as white paper under light conditions similar to "; brightness(x);

cout << endl << endl << endl; }

void testData(char* data,double CONST, char* output) { Function a; a.getData(data); a.findLambda(CONST);

ofstream fout(output);

if(!fout) {cout << "Could not open file " << output << " program terminated." << endl; fout.close(); return;}

for(double i=MINI; i<=MAXI ; i+=0.1) fout << i << "\t" << a.value(i) << endl;

// for(double i=0; i<a.size() ; i++) // cout << a.wat(i) << "\t" << a.fat(i)<<"\t" << a.lat(i)<< endl;

fout.close();

}

double luminousFlux(Function x, Function y, double cden, double turnover) { double sum1=0, sum2=0; int a=350,b=750; double h=H(a,b,N);

for(double j=1; j<=(N/2)-1;j++) {sum1+=x.value(X(2*j,a,h))*y.value(X(2*j,a,h))/X(2*j,a,h);}

for(double j=1; j<=N/2; j++) {sum2+=x.value(X(2*j-1,a,h))*y.value(X(2*j-1,a,h))/X(2*j-1,a,h);}

return 683.002/(4*3.141)*HCON*CCON*turnover*cden*1e9*((h/3)* ( x.value(a)*y.value(a)/a+ 2*sum1 + 4*sum2 + x.value(b)*y.value(b)/b ))/integral(a,b,y);

}

double integral (double a, double b, Function func) { double sum1=0, sum2=0; double h=H(a,b,N);

for(double j=1; j<=(N/2)-1;j++) {sum1+=func.value( X( 2*j , a , h ) );}

for(double j=1; j<=N/2; j++) {sum2+=func.value( X( (2*j)-1 , a , h ) );}

return (h/3)* ( func.value(a)+ 2*sum1 + 4*sum2 + func.value(b) ); }

//.............................................................................................................

void Function::getData(const char* a) { double num;

SIZE=200;

ifstream fin(a); if(!fin) {cout << "Could not open file " << a << " program terminated." << endl; fin.close(); return;}

fin >> num; for (int i=0; !fin.eof(); i++ ) { if(!(i%2)){ wavelength.push_back(num); fin >> num;} else { function.push_back(num); fin >> num;} }

fin.close();

int total=wavelength.size(); vector <double> tempwavelength; vector <double> tempfunction;

if(total>SIZE){ for(int i=0; i<wavelength.size(); i++) { if(i%(total/SIZE)==0){ tempwavelength.push_back(wavelength.at(i)); tempfunction.push_back(function.at(i)); }} wavelength.clear(); function.clear(); for(int i=0; i<tempwavelength.size(); i++) { wavelength.push_back(tempwavelength.at(i)); function.push_back(tempfunction.at(i)); } }

width=(wavelength.back()-wavelength.front())/wavelength.size();

}

void Function::print(double min, double max, double increment) { for(double i=min; i<=max ; i+=increment) cout << i << "\t" << value(i) << endl; }

void Function::printStore() { for(double i=0; i<wavelength.size() ; i++) cout << wavelength.at(i) << "\t" << function.at(i) << endl;

}

void Function::print() { for(int i=0; i<wavelength.size(); i++) { cout << wavelength.at(i) << "\t" << function.at(i) << endl; }

}

double Function::RBF(double a, double b) { return exp(-(a-b)*(a-b)/(2*sigma*sigma)); //return fabs(a-b); //return 1/sqrt(1+sigma*(a-b)*(a-b)); }

void Function::findLambda(double CONST) { sigma=width*CONST;

int M=wavelength.size();

 matrix <double> PHI(M,M);
 matrix <double> lamb(M,1);
 matrix <double> func(M,1);  

for (int i=0; i<M; i++) { func.setvalue(i,1,function.at(i)); for (int j=0; j<M; j++) { PHI.setvalue(i,j,RBF(wavelength.at(i),wavelength.at(j)));

}}

PHI.invert(); lamb.settoproduct(PHI,func);

double value;bool success;

for(int i=0; i<M; i++){lambda.push_back(lamb.getvalue(i,1));} }

double Function::value(double x) { double sum=0; for(int i=0; i<wavelength.size(); i++) sum+=lambda.at(i)*RBF(x,wavelength.at(i));

return sum; }

void Function::fix() { for(int i=0; i<wavelength.size(); i++){ for(int j=0;j<i; j++){ if(wavelength.at(i)==wavelength.at(j)) {wavelength.erase(wavelength.begin()+j); function.erase(function.begin()+j); j=j-1;} }}

}