Algorithms

from numpy import *
import sys
import numpy as np
x=1
epsilon=sys.float_info.epsilon/2.0
# machine epsilon
h = np.linspace(1,epsilon,1000)
e1= (((sin(x)- sin(x-h))/h-cos(x))/cos(x))                      #realtive error in first order
e2= (((sin(x-2*h)-4*sin(x-h)+3*sin(x))/(2*h)-cos(x))/cos(x))    #realtive error in second order
e3= (((2*sin(x+h)+3*sin(x)-6*sin(x-h)+sin(x-2*h))/(6*h)-cos(x))/cos(x))#realtive error in third order
e4 = (((-sin(x-3*h)+6*sin(x-2*h)-18*sin(x-h)+10*sin(x)+3*sin(x+h))/(12*h)-cos(x))/cos(x))#realtive error in fourth order
print(h)
import matplotlib.pyplot as plt
plt.plot(h,e1,'-b',lw=3)
plt.plot(h,e2,'-g',lw=3)
plt.plot(h,e3,'-m',lw=3)
plt.plot(h,e4,'-r',lw=3)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('h')
plt.ylabel('Relative error')
plt.legend(['first order','second order','third order','fourth order'])
plt.title('Effect of order of accuracy on the numerical differentiation')
plt.show()

template <typename T>
constexpr MCSRmatrix<T>::MCSRmatrix( std::initializer_list<std::initializer_list<T>> rows)
{
      this->dim  = rows.size();
      auto _rows = *(rows.begin());
            
      aa_.resize(dim+1);
      ja_.resize(dim+1);
      
      if(dim != _rows.size())
      {
          throw InvalidSizeException("Error in costructor! MCSR format require square matrix!");  
      }
      
      itype w = 0 ;
      ja_.at(w) = dim+2 ;
      for(auto ii = rows.begin(), i=1; ii != rows.end() ; ++ii, i++)
      {
          for(auto ij = ii->begin(), j=1, elemCount = 0 ; ij != ii->end() ; ++ij, j++ )   
          {
              if(i==j)
                 aa_[i-1] = *ij ;
              else if( i != j && *ij != 0 )
              {   
                 ja_.push_back(j); 
                 aa_.push_back(*ij); 
                 elemCount++ ;
              }
              ja_[i] = ja_[i-1] + elemCount;           
          }
      }     
}

template <typename T>
std::size_t constexpr MCSRmatrix<T>::findIndex(const itype row ,  const itype col) const noexcept 
{    
     assert( row > 0 && row <= dim && col > 0 && col <= dim ); 
     if(row == col)
     {
       return row-1;
     }
     int i = -1;
     
      for(int i = ja_.at(row-1)-1 ; i < ja_.at(row)-1 ; i++ )
      {
            if( ja_.at(i) == col ) 
            { 
                  
              return i   ;
            }
      }
      return -1;
}

template <typename T>
const T MCSRmatrix<T>::operator()(const itype r , const itype c) const noexcept 
{     
     auto i = findIndex(r,c); 
     if( i != -1 && i < aa_.size() )
     {      
            return aa_.at(i) ;
     }
     else
     {
            return 0.0 ;
     }       
}

template<typename T>
auto constexpr MCSRmatrix<T>::insertAt(itype row ,itype col , T elem) noexcept 
{
      if(elem == 0)
      {
          std::cerr << "zero element to insert ! Exit 1" << std::endl; 
      }
      
      int index = findIndex(row,col);
      if( index >= dim+1  || row == col)
      {     
          aa_.at(index) = elem ;
      }
      else if(index == -1) 
      {

          for(auto i=row; i< dim+1 ; i++) 
          {     
                ja_.at(i) += 1 ;
          }
          
          if(ja_.at(row) >= aa_.size() ) //std::cout << aa_.size() << ' ' <<ja_.at(row) << std::endl;
          {
              aa_.push_back(elem);
              ja_.push_back(col);
          }
          else if ( ja_.at(row) < aa_.size() )
          {
                ja_.insert(ja_.begin() + dim  + row , col );
                aa_.insert(aa_.begin() + dim  + row , elem);
          }
      }
}

// interval = milliseconds, seconds, minutes, hours, days, months, year
//amount (interval) is the amount you need to add to the timestamp in a certain interval.  (eg amount = 12, interval = seconds ==> add 12 seconds)
function FindGaps(DateTime dts_start, DateTime dts_end, int amount, string interval, Array list2check){

//Build reference list.
Array listreference<DateTime, boolean>;     //dictionary containing datetime as key and boolean as value (or integer, enums, ...)
for(DateTime dts = dts_start; dts < dts_end; dts += amount (interval)){
    listreference.Add(dts, false);      //initialize base values
}

//check list with reference list.
for(int i = 0; i < list2check.length; i++){     //or use long ?
    //if value of DateTime list2check is in the reference list:  Jeej!
    if(KeyExists(listrefence, list2check[i])){
        listrefence[list2check[i]]) = true;
    }
}

return listrefence;
}