@lcnature, would you like me to go over it now?

@CameronTEllis That would be awesome. One change I made was in fmrisim.py, which especially needs your eyes on it in case I make new bugs.
That change in fmrisim.py was primarily due to some strange corner case I encountered, in which machine precision in float number cause that conditional clause to execute when it is actually unnecessary (eg., sometimes total_time * temporal_resolution can be smaller than its integer version by a very small amount due to machine precision, and then int(total_time * temporal_resolution) somehow is smaller than int(round(total_time * temporal_resolution)) by 1 in very occasional case...).

And then I realized that the conditional clause in generate_stimfunction to shorten the stimfunction may not be necessary, because the lines above won't actually lengthen it as Matlab would do even if offset_idx happens to be larger than int(round(total_time * temporal_resolution)).
If you want to keep that clause, Line 528 might need to be changed to stimfunction = stimfunction[0:int(total_time * temporal_resolution), :], otherwise stimfunction can become 1-d nparray instead of a 2-d one, which convolve_hrf is not happy about.