Making a fitfunction multithreaded

Average rating
(1 vote)

Sometimes you will want to make your fit function multithreaded. This is not much use when you use Funcfit or Curvefit, as the operation does the threading itself. However, it can be very useful when using GenCurvefit as this XOP operation cannot perform multithreading itself, so you need to do it in the fitfunction. The code list below is a demonstration of how to do this. Simply change the code in the _kernel functions and you're off. You can of course change all the function names to your hearts content.

You should note that significant speedups are only obtained when the fitfunction to be calculated takes a long time. If it is too quick then the thread overhead will kill you. Test the speedups before you use code exclusively

Function myfitfunc_wrapper(w, yy, xx) : fitfunc
Wave w, yy, xx	
	Multithread yy = myfitfunc_kernel(w, xx)
End
 
Threadsafe function myfitfunc_kernel(w, xx)
Wave w
variable xx
	//your code goes here
	return w[0] + sqrt(xx) * w[1] * sin(xx) * cos(x) * ln(x)
End

-OR-
Function myfitfunc_AAO_wrapper(w, yy, xx) : fitfunc
Wave w, yy, xx
 
	variable nt = ThreadProcessorCount
	variable tgID = ThreadGroupCreate(nt)
	variable ii, startP, perT, remainingP
 
	perT = floor(dimsize(yy, 0) / nt)
	remainingP = dimsize(yy, 0)
 
	make/n=(nt)/free/wave theDataY, theDataX
	for(ii = 0 ; ii < nt ; ii+=1)
		if(ii < nt -1)
			theDataY[ii] = newfreewave(0x04, perT)
			theDataX[ii] = newfreewave(0x04, perT)
		else
			theDataY[ii] = newfreewave(0x04, remainingP)
			theDataX[ii] = newfreewave(0x04, remainingP)		
		endif
		Wave yyy = theDataY[ii]
		Wave xxx = theDataX[ii]
 
		yyy = yy[p+startP]
                xxx = xx[p+startP]
 
		startP += perT
		remainingP -= perT
		ThreadStart tgID, ii, myfitfunc_AAO_kernel(w, theDataY[ii], theDataX[ii])
		WaveClear yyy, xxx
	endfor
	variable finished = ThreadGroupWait(tgID, inf)
	finished = threadgrouprelease(tgID)
 
	//now copy back into structure wave	
	startP = 0
	perT = floor(dimsize(yy, 0) / nt)
	remainingP = dimsize(yy, 0)
	for(ii = 0 ; ii < nt ; ii+=1)
		Wave yyy = theDataY[ii]
		yy[startP, startP+numpnts(yyy)] = yyy[p-startP]
		Waveclear yyy
 
		startP += perT
		remainingP -= perT
	endfor
End
 
Threadsafe function myfitfunc_AAO_kernel(w, yy, xx)
Wave w, yy, xx
	//your code goes here
	yy = w[0] + w[1] * sqrt(xx) * sin(xx) * cos(x) * ln(x)
End

-OR- if you wanted to do a structure fit

Structure fitfuncStruct   
	Wave w
	wave y
	wave x[50]
 
	int16 numVarMD
	wave ffsWaves[50]
	wave ffsTextWaves[10]
	variable ffsvar[5]
	string ffsstr[5]
	nvar ffsnvars[5]
	svar ffssvars[5]
	funcref allatoncefitfunction ffsfuncrefs[10]
	uint32 ffsversion    // Structure version. 
	EndStructure 
 
Threadsafe Function allatoncefitfunction(w,y,x)
	Wave w,y,x
	//you need this to act as a signature for the fitfuncStruct
	//it needs to be threadsafe
End
 
Threadsafe function myfitfunc_AAO_kernel(w, yy, xx)
Wave w, yy, xx
	//your code goes here
	yy = w[0] + w[1] * sqrt(xx) * sin(xx) * cos(x) * ln(x)
End
 
Function myfitfunc_AAO_wrapper(s) : fitfunc
	Struct fitfuncStruct &s
 
 	//work out how many points per thread
	variable nt =  ThreadProcessorCount
	variable tgID = ThreadGroupCreate(nt)
	variable ii, startP, perT, remainingP, finished
	perT = floor(dimsize(s.y, 0) / nt)
	remainingP = dimsize(s.y, 0)
 
 	//create waveref waves to hold references to the y and x waves to be distributed amongst each thread.
 	//you could create more than y and x if your function took more arguments
	make/n=(nt)/free/wave theDataY, theDataX
 
	//create the coefficients wave to be sent to all the threads, as you can't pass in s.w to a thread
	Wave theCoefs = s.w
 
	//this is the function we are going to pass the coefficients, and y and x waves onto. 
	 //-It doesn't have to be an allatoncefitfunction it can be any function type you want (so long as it doesn't take structures).
	//- It needs to be threadsafe!
	//- the function can take as many parameters as you want, but you'll have to create a free wave for each of them (like theCoefs) to 
	//    be passed into the thread.
	Funcref allatoncefitfunction theFitFunction = s.ffsfuncrefs[0]
 
	for(ii = 0 ; ii < nt ; ii+=1)
		if(ii < nt -1)
			theDataY[ii] = newfreewave(0x04, perT)
			theDataX[ii] = newfreewave(0x04, perT)
		else
			theDataY[ii] = newfreewave(0x04, remainingP)
			theDataX[ii] = newfreewave(0x04, remainingP)		
		endif
		Wave yyy = theDataY[ii]
		Wave xxx = theDataX[ii]
 
		yyy = s.y[p+startP]
		xxx = s.x[0][p+startP]
 
		startP += perT
		remainingP -= perT
 
		//NOTE that you can't pass in a structure parameter in a threadstart operation. See:
		//DisplayHelptopic "Threadstart"
		// for more information
		ThreadStart tgID, ii, myfitfunc_AAO_kernel(theCoefs, yyy, xxx)
		WaveClear yyy, xxx
	endfor
	 finished = ThreadGroupWait(tgID, inf)
	finished = threadgrouprelease(tgID)
 
	//now copy back into structure wave	
	startP = 0
	perT = floor(dimsize(s.y, 0) / nt)
	remainingP = dimsize(s.y, 0)
	for(ii = 0 ; ii < nt ; ii+=1)
		Wave yyy = theDataY[ii]
		s.y[startP, startP+numpnts(yyy)] = yyy[p-startP]
		Waveclear yyy
 
		startP += perT
		remainingP -= perT
	endfor
End

Hi, first of all thanks for

Hi, first of all thanks for your post! I am using pretty massive all-at-once fitting functions, and I wanted to see how they perform using genetic optimisation code. I have an i7 CPU, so I want to exploit multithreading. I tried your code but it didn't work, as I received the error "FitFunction returned NaN/INF for one or more values". I specify that these fitfunctions in their "normal" form are working flawlessly. Do you know the possible source of error? I attach a .txt of my code. Let me know if you need further info.
Thanks!

Riccardo

AttachmentSize
fitfunc_code.txt4.19 KB

can you post an experiment

can you post an experiment with some data?

Can't promise I can solve your problem, but I'll have a look.

I'll shrink the one I am

I'll shrink the one I am using right now, to make it lighter. As soon as I can I will post it!

Back to top