source: issm/oecreview/Archive/13393-13976/ISSM-13568-13569.diff@ 27230

../trunk-jpl/src/c/shared/Sorting/binary_search.cpp

@@ -10 +10 @@
 
 #include <stdio.h>
 
-int binary_search(int* poffset,int target,int* sorted_integers,int num_integers){
+int binary_search(int* poffset,int target,int* sorted_integers,int num_integers){ /*{{{*/
 
         /*output: */
         int offset;  //offset, if found
@@ -62 +62 @@
         
         /*Return result: */
         return found;
-}
+} /*}}}*/
+int binary_search(int* poffset,double target,double* sorted_doubles,int num_doubles){ /*{{{*/
+
+        /*output: */
+        int offset=0;  
+        int found=0; /*found=0: not found.   
+                                   found=1: found, and target is == to value at offset 
+                                   found=2: found, and target is > to value at offset  and < to value at offset+1
+                                   */
+
+        /*intermediary: */
+        double* beg=NULL;
+        double* end=NULL;
+        double* mid=NULL;
+
+        // point to beginning and end of the array
+        beg=sorted_doubles;
+        end=sorted_doubles+num_doubles;
+        mid=beg+(int)(num_doubles/2.0);
+
+        if (target<*beg){
+                offset==-1;
+                found=0;
+        }
+        if (*beg==target){
+                found=1;
+                offset=0;
+        }
+        else if(*(end-1)==target){
+                found=1;
+                offset=num_doubles-1;
+        }
+        else{
+                while((beg <= end) && !( target>=*mid &&  target<*(mid+1)) ){
+                        // is the target in lower or upper half?
+                        if (target < *mid) {
+                                end = mid - 1;  //new end
+                                mid = beg + (end-beg)/2;  //new middle
+                        }
+                        else {
+                                beg = mid + 1;  //new beginning
+                                mid = beg + (end-beg)/2;  //new middle
+                        }
+                }
+                          
+                //did we find the target?
+                if( target>*mid &&  target<*(mid+1)){
+                        found=2;
+                        offset=mid-sorted_doubles;
+                }
+                else if( target==*mid){
+                        found=1;
+                        offset=mid-sorted_doubles;
+                }
+                else {
+                        found=0;
+                }
+        }
+
+        /*Assign output pointers:*/
+        *poffset=offset;
+        
+        /*Return result: */
+        return found;
+} /*}}}*/

../trunk-jpl/src/c/shared/Sorting/sorting.h

@@ -6 +6 @@
 #define  _SORTING_H_
 
 int binary_search(int* poffset,int target,int* sorted_integers,int num_integers);
+int binary_search(int* poffset,double target,double* sorted_doubles,int num_doubles);
 
 #endif //ifndef _SORTING_H_

../trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp

@@ -425 +425 @@
 /*FUNCTION TransientInput::GetTimeInput{{{*/
 Input* TransientInput::GetTimeInput(IssmDouble intime){
 
-        int     i,j;
         IssmDouble  deltat;
         IssmDouble  alpha1,alpha2;
-        bool    found=false;
+        
         Input*  input=NULL;
         Input*  input1=NULL;
         Input*  input2=NULL;
+        
+        bool    found=false;
+        int     found_offset=0;
+        int     offset;
 
         /*Ok, we have the time, go through the timesteps, and figure out which interval we 
          *fall within. Then interpolate the values on this interval: */
-        if(intime<this->timesteps[0]){
-                /*get values for the first time: */
-                input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
-                found=true;
+        found_offset=binary_search(&offset,intime,this->timesteps,this->numtimesteps);
+
+        if (found_offset==0){
+                if(intime<this->timesteps[0]){
+                        /*get values for the first time: */
+                        input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
+                        found=true;
+                }
+                else if(intime>this->timesteps[this->numtimesteps-1]){
+                        /*get values for the last time: */
+                        input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy();
+                        found=true;
+                }
         }
-        else if(intime>this->timesteps[this->numtimesteps-1]){
-                /*get values for the last time: */
-                input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy();
-                found=true;
-        }
         else{
-                /*Find which interval we fall within: */
-                for(i=0;i<this->numtimesteps;i++){
-                        if(intime==this->timesteps[i]){
-                                /*We are right on one step time: */
-                                input=(Input*)((Input*)this->inputs->GetObjectByOffset(i))->copy();
-                                found=true;
-                                break; //we are done with the time interpolation.
-                        }
-                        else{
-                                if(this->timesteps[i]<intime && intime<this->timesteps[i+1]){
-                                        /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
-                                        deltat=this->timesteps[i+1]-this->timesteps[i];
-                                        alpha2=(intime-this->timesteps[i])/deltat;
-                                        alpha1=(1.0-alpha2);
+                if (found_offset==1){
+                        /*intime is exactly equal to this->timesteps[offset]: */
+                        input=(Input*)((Input*)this->inputs->GetObjectByOffset(offset))->copy();
+                        found=true;
+                }
+                else if (found_offset==2){
+                        
+                        /*ok, we have the interval ]offset:offset+1[. Interpolate linearly for now: */
+                        deltat=this->timesteps[offset+1]-this->timesteps[offset];
+                        alpha2=(intime-this->timesteps[offset])/deltat;
+                        alpha1=(1.0-alpha2);
 
-                                        input1=(Input*)this->inputs->GetObjectByOffset(i); 
-                                        input2=(Input*)this->inputs->GetObjectByOffset(i+1);
+                        input1=(Input*)this->inputs->GetObjectByOffset(offset); 
+                        input2=(Input*)this->inputs->GetObjectByOffset(offset+1);
 
-                                        input=(Input*)input1->copy();
-                                        input->Scale(alpha1);
-                                        input->AXPY(input2,alpha2);
+                        input=(Input*)input1->copy();
+                        input->Scale(alpha1);
+                        input->AXPY(input2,alpha2);
 
-                                        found=true;
-                                        break;
-                                }
-                                else continue; //keep looking on the next interval
-                        }
+                        found=true;
                 }
+                else _error_("binary_search returned the following value for found: " << found_offset<< "which is not supported yet!");
         }
+
         if(!found)_error_("did not find time interval on which to interpolate forcing values!");
 
         /*Assign output pointer*/