A custom fitting equation that also calls a fit

SailBlue5 wrote: I am using version 6.2.2.2, downloaded recently. It doesn't give any warning messages. I feared that this was the case.

Any work arounds I could do?

#pragma rtGlobals=1	// Use modern global access method.

//This function works as a basic fit function for Diffusion's dependence on oscillating amplitude 
Function DiffusionVsOscAmplitude(w,OscAmp) : FitFunc
	Wave w
	Variable OscAmp

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ This runs a simulation to the differential equation determine the fit.
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ OscAmp
	//CurveFitDialog/ Coefficients 2
	//CurveFitDialog/ w[0] = Ds
	//CurveFitDialog/ w[1] = invtm
print "Starting fitfunc"

	//Set global variables
	variable /G position_cm = 200e-7
	variable /G SimLength_cm =400e-7
	Variable /G dx_cm = 2e-7
	variable /G dt_s = .0002
	variable /G SimTime_s = 3
	
	//Make the sensitive volume
	SenSliceFull(OscAmp)
	print OscAmp
	print w[0]
	print w[1]
print "Calling runSim"	
	return RunSimulation(w[0],w[1])
End


//To run the simulation 
Function RunSimulation(Ds, invtm)
	Variable Ds     //diffusion constant
	Variable invtm  //Relaxation rate for T_1 or \tau_m
print "In RunSim"
	variable /G position_cm
	variable /G SimLength_cm
	Variable /G dx_cm
	variable /G dt_s
	variable /G SimTime_s
	Wave ExpFits
	Wave SenSlice
	Wave Spins
	
	variable Iterations = round(SimTime_s/dt_s) //Number of spatial points to simulate
	variable Positions = round(SimLength_cm/dx_cm)  //Number of time steps to take
	variable frac = cosh(dt_s*invtm)*exp(-1*dt_s*invtm)  //the fraction of spins that flip in time dt.
	variable i  //Will hold the time step
	variable j //will hold the position step
	variable k //used to update plots every now and then
	variable sum  //will hold the integration of Spin*SenSlice
	variable result  //The resultant measured inverse corrolation time of the force
	
	Make /O /N=(Positions) TempSpins  //This will hold the spin state at the next time step temporarly
	Make /O /N=(Iterations) OrigSpinState  //This will hold how much of the original spin state is left convolved with the sen slice vs time (for each dt)
	SetScale/P x 0,(dt_s),"", OrigSpinState  //each point repersents a single time step
	Make /O /N=2 TempCoef  //Since we will have a fit inside a fit we need to do house cleaning so the two W-coef waves do not screw with each other
	
	OrigSpinState = nan  //just so I don't have lingering values from a previous run
	for(i=0;i<Iterations;i=i+1)
		sum = 0
		for(j=0;j<Positions; j=j+1)
			sum = sum + SenSlice[j]*Spins[j]  //Integrate to find "force"
			
			//Spins diffusing
			if(j==0)
				TempSpins[j] = Spins[j] + Ds*dt_s*(Spins[j+1]-Spins[j])/(dx_cm*dx_cm)  	//spins are blocked from diffusing past on the left side
				//TempSpins[j] = 0													//Spins diffuse past the left side (and disapear)	
			elseif(j==Positions-1)
				TempSpins[j] = Spins[j] + Ds*dt_s*(Spins[j-1]-Spins[j])/(dx_cm*dx_cm)  //spins are blocked from diffusing past on the right side
				//TempSpins[j] = 0													//Spins diffuse past the right side (and disapear)
			else
				TempSpins[j] = Spins[j] + Ds*dt_s*(Spins[j-1]+Spins[j+1]-2*Spins[j])/(dx_cm*dx_cm)  
			endif
		endfor
		
		OrigSpinState[i] = sum  //The "force" at this time period
		
		for(j=0;j<Positions; j=j+1)
			Spins[j] = frac*TempSpins[j]  //Get ready for next iteration
		endfor
		if(k<1)
			doUpdate
			k=Iterations/100
		endif
		k=k-1
	endfor
	print "Calling fit"
	CurveFit/NTHR=0 exp kwCWave=ExpFits,  OrigSpinState /D 
	print ExpFits[2]
	print "Returning"
	return ExpFits[2]
end

Function SenSliceFull(OscAmp)
	variable OscAmp
	variable /G position_cm
	variable /G SimLength_cm
	Variable /G dx_cm
	
	variable i
	
	Make /O /N=(round(SimLength_cm/dx_cm)) SenSlice
	Make /O /N=(round(SimLength_cm/dx_cm)) Spins
	for(i=0;i<round(SimLength_cm/dx_cm);i=i+1)
		SenSlice[i]=0
		Spins[i]=0
	endfor
	for(i=ceil(-1*(OscAmp/dx_cm));i<=floor(OscAmp/dx_cm);i=i+1)
		SenSlice[round(position_cm/dx_cm)+i]=Cos((i*dx_cm/OscAmp)*(pi/2))
		Spins[round(position_cm/dx_cm)+i]=1
	endfor
	SetScale/P x 0,(dx_cm),"", SenSlice
	SetScale/P x 0,(dx_cm),"", Spins
end

threadsafe Function DoMyFit(cwave, datawave)
	Wave cwave, datawave
	
	CurveFit/NTHR=1 exp kwCWave=cwave,  datawave
end

Function runsimulation(Ds, invtm)
	Variable Ds     //diffusion constant
	Variable invtm  //Relaxation rate for T_1 or \tau_m
	
	//*** lots of code here****
	
	Variable threadid = ThreadGroupCreate(1)		// just one thread to encapsulate the inner fit
	ThreadStart threadid, 0, DoMyFit(ExpFits,  OrigSpinState)
	do
		variable tgs= ThreadGroupWait(threadid,5000)		// five-second wait is excessive, but we don't want to spin the loop.
												// We just want to wait until the inner fit is done, which could be fairly lengthy.
	while( tgs != 0 )

	return ExpFits[2]
end

johnweeks wrote: Michael-

Here is some code that might serve as a start on trying my idea. To read about threading, and what the calls in this code do, execute this command:

DisplayHelpTopic "ThreadSafe Functions and Multitasking"

Note that my code is not complete, so I haven't tested it. But it should serve as a starting point. I will be very interested to know if this works!

threadsafe Function DoMyFit(cwave, datawave)
	Wave cwave, datawave
	
	CurveFit/NTHR=1 exp kwCWave=cwave,  datawave
end

Function runsimulation(Ds, invtm)
	Variable Ds     //diffusion constant
	Variable invtm  //Relaxation rate for T_1 or \tau_m
	
	//*** lots of code here****
	
	Variable threadid = ThreadGroupCreate(1)		// just one thread to encapsulate the inner fit
	ThreadStart threadid, 0, DoMyFit(ExpFits,  OrigSpinState)
	do
		variable tgs= ThreadGroupWait(threadid,5000)		// five-second wait is excessive, but we don't want to spin the loop.
												// We just want to wait until the inner fit is done, which could be fairly lengthy.
	while( tgs != 0 )

	return ExpFits[2]
end

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

SailBlue5 wrote: Creating a new thread worked. I was able to perform the fitting function within a fitting function.