#!/usr/bin/env python
#
#   Copyright 2016 Blaise Frederick
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.
#
#
#       $Author: frederic $
#       $Date: 2016/06/14 12:04:51 $
#       $Id: simdata,v 1.18 2016/06/14 12:04:51 frederic Exp $
#
from __future__ import print_function
import sys
import rapidtide.miscmath as tide_math
import rapidtide.io as tide_io
import rapidtide.resample as tide_resample
from matplotlib.pyplot import *


def usage():
    print(
        "usage: simdata fmrifilename immeanfilename boldpcfilename lagfilename regressor regressorsamprate numskip inputstarttime outputname noiselevel sliceorder")
    print("")
    print("required arguments:")
    print("	fmrifilename		- the example BOLD fmri file")
    print("	immeanfilename		- the 3D mean image file")
    print("	boldpcfilename		- the 3D bold percentage file")
    print("	lagfilename		- the 3D time delay (in seconds) at every point")
    print("	regressorfile		- the name of the regressor text file")
    print("	regressorsamprate	- the sample rate of the regressor file, in Hz")
    print("	numskip			- the use to simulate tr periods deleted during preprocessing")
    print(
        "	inputstarttime		- the time delay, in seconds, into the inputfile that matches the start time of the fmrifile")
    print("	outputname		- the root name for the output files")
    print("	noiselevel		- the variance of the noise")
    print("	sliceorder		- slice acquisition order:")
    print("					0 : None")
    print("					1 : Regular up (0, 1, 2, 3, ...)")
    print("					2 : Regular down")
    print("					3 : Use slice order file")
    print("					4 : Use slice timings file")
    print("					5 : Standard Interleaved (0, 2, 4 ... 1, 3, 5 ... )")
    print("					6 : Siemens Interleaved (0, 2, 4 ... 1, 3, 5 ... for odd number of slices)")
    print("					                        (1, 3, 5 ... 0, 2, 4 ... for even number of slices)")

    return ()


# set default variable values
displayplots = False
fastresample = True
debug = False

# get the command line parameters
if len(sys.argv) != 12:
    usage()
    exit()

# handle required args first
fmrifilename = sys.argv[1]
immeanfilename = sys.argv[2]
boldpcfilename = sys.argv[3]
lagfilename = sys.argv[4]
filename = sys.argv[5]
inputfreq = float(sys.argv[6])
fmriskip = int(sys.argv[7])
inputstarttime = float(sys.argv[8])
outputname = sys.argv[9]
noiselevel = float(sys.argv[10])
sliceordertype = int(sys.argv[11])

print("fmrifilename", fmrifilename)
print("immeanfilename", immeanfilename)
print("boldpcfilename", boldpcfilename)
print("lagfilename", lagfilename)
print("filename", filename)
print("inputfreq", inputfreq)
print("fmriskip", fmriskip)
print("inputstarttime", inputstarttime)
print("outputname", outputname)
print("noiselevel", noiselevel)

fmritr, numtrs = tide_io.fmritimeinfo(fmrifilename)
nim_fmri, fmridata, fmriheader, fmridims, fmrisizes = tide_io.readfromnifti(fmrifilename)
print("fmri data: ", numtrs, " timepoints, tr = ", fmritr)

# read in the timecourse to resample
inputvec = tide_math.stdnormalize(tide_io.readvec(filename))
simregressorpts = len(inputvec)

# prepare the input data for interpolation
print("Input regressor has ", simregressorpts, " points")
inputstep = 1.0 / inputfreq
nirs_x = r_[0.0:1.0 * simregressorpts] * inputstep - inputstarttime
nirs_y = inputvec[0:simregressorpts]
print('nirs regressor runs from ', nirs_x[0], ' to ', nirs_x[-1])

# prepare the output timepoints
fmrifreq = 1.0 / fmritr
initial_fmri_x = r_[0.0:fmritr * (numtrs - fmriskip):fmritr] + fmritr * fmriskip
print('length of fmri after removing skip:', len(initial_fmri_x))
print('fmri time runs from ', initial_fmri_x[0], ' to ', initial_fmri_x[-1])

# read in the immean, boldpc and lag files
print("reading in source files")
nim_immean, immeandata, immeanheader, immeandims, immeansizes = tide_io.readfromnifti(immeanfilename)
if not tide_io.checkspacedimmatch(immeandims, fmridims):
    print("immean file does not match")
    exit()
nim_boldpc, boldpcdata, boldpcheader, boldpcdims, boldpcsizes = tide_io.readfromnifti(boldpcfilename)
if not tide_io.checkspacedimmatch(boldpcdims, fmridims):
    print("boldpc file does not match")
    exit()
nim_lag, lagdata, lagheader, lagdims, lagsizes = tide_io.readfromnifti(lagfilename)
if not tide_io.checkspacedimmatch(lagdims, fmridims):
    print("lag file does not match")
    exit()

# now set up the simulated data array
thesizes = fmrisizes
xdim = thesizes[0]
ydim = thesizes[1]
slicethickness = thesizes[2]
tr = thesizes[3]

thedims = fmridims
xsize = thedims[1]
ysize = thedims[2]
numslices = thedims[3]
timepoints = thedims[4]

simdata = zeros((xsize, ysize, numslices, len(initial_fmri_x)), dtype='float')

fmrilcut = 0.0
fmriucut = fmrifreq / 2.0

# set up fast resampling
padvalue = 60.0
numpadtrs = int(padvalue / fmritr)
padvalue = fmritr * numpadtrs

genlagtc = tide_resample.fastresampler(nirs_x, nirs_y, padvalue=padvalue, doplot=False)
fmri_y_f = genlagtc.yfromx(initial_fmri_x)
fmri_y = tide_resample.doresample(nirs_x, nirs_y, initial_fmri_x)
# fig = figure()
# ax = fig.add_subplot(111)
# plot(initial_fmri_x, fmri_y_f + 1.0, initial_fmri_x, fmri_y)
# show()

# fig = figure()
# ax = fig.add_subplot(111)
# ax.set_title('Regressors')
# plot(nirs_x, nirs_y)
# show()

# loop over space
for k in range(0, numslices):
    sliceoffsettime = tide_resample.calcsliceoffset(sliceordertype, k, numslices, fmritr)
    print("processing slice ", k, ": sliceoffsettime=", sliceoffsettime)
    for j in range(0, ysize):
        for i in range(0, xsize):
            # correct the output time points
            fmri_x = initial_fmri_x - lagdata[i, j, k] - sliceoffsettime
            if debug:
                print(i, j, k, fmri_x[0], initial_fmri_x[0], lagdata[i, j, k], sliceoffsettime)

            if fastresample:
                fmri_y = genlagtc.yfromx(fmri_x)
            else:
                fmri_y = tide_resample.doresample(nirs_x, nirs_y, fmri_x)
            thenoise = noiselevel * standard_normal(len(fmri_y))
            simdata[i, j, k, :] = immeandata[i, j, k] * (1.0 + (boldpcdata[i, j, k] / 100.0) * fmri_y) + thenoise

tide_io.savetonifti(simdata, fmriheader, fmrisizes, outputname)

994d5aba774dc5c0a7c9e2f5ce0180aabd172a5b
<mask token>

<mask token>


def usage():
    print(
        'usage: simdata fmrifilename immeanfilename boldpcfilename lagfilename regressor regressorsamprate numskip inputstarttime outputname noiselevel sliceorder'
        )
    print('')
    print('required arguments:')
    print('\tfmrifilename\t\t- the example BOLD fmri file')
    print('\timmeanfilename\t\t- the 3D mean image file')
    print('\tboldpcfilename\t\t- the 3D bold percentage file')
    print('\tlagfilename\t\t- the 3D time delay (in seconds) at every point')
    print('\tregressorfile\t\t- the name of the regressor text file')
    print('\tregressorsamprate\t- the sample rate of the regressor file, in Hz'
        )
    print(
        '\tnumskip\t\t\t- the use to simulate tr periods deleted during preprocessing'
        )
    print(
        '\tinputstarttime\t\t- the time delay, in seconds, into the inputfile that matches the start time of the fmrifile'
        )
    print('\toutputname\t\t- the root name for the output files')
    print('\tnoiselevel\t\t- the variance of the noise')
    print('\tsliceorder\t\t- slice acquisition order:')
    print('\t\t\t\t\t0 : None')
    print('\t\t\t\t\t1 : Regular up (0, 1, 2, 3, ...)')
    print('\t\t\t\t\t2 : Regular down')
    print('\t\t\t\t\t3 : Use slice order file')
    print('\t\t\t\t\t4 : Use slice timings file')
    print('\t\t\t\t\t5 : Standard Interleaved (0, 2, 4 ... 1, 3, 5 ... )')
    print(
        '\t\t\t\t\t6 : Siemens Interleaved (0, 2, 4 ... 1, 3, 5 ... for odd number of slices)'
        )
    print(
        '\t\t\t\t\t                        (1, 3, 5 ... 0, 2, 4 ... for even number of slices)'
        )
    return ()


<mask token>
if len(sys.argv) != 12:
    usage()
    exit()
<mask token>
print('fmrifilename', fmrifilename)
print('immeanfilename', immeanfilename)
print('boldpcfilename', boldpcfilename)
print('lagfilename', lagfilename)
print('filename', filename)
print('inputfreq', inputfreq)
print('fmriskip', fmriskip)
print('inputstarttime', inputstarttime)
print('outputname', outputname)
print('noiselevel', noiselevel)
<mask token>
print('fmri data: ', numtrs, ' timepoints, tr = ', fmritr)
<mask token>
print('Input regressor has ', simregressorpts, ' points')
<mask token>
print('nirs regressor runs from ', nirs_x[0], ' to ', nirs_x[-1])
<mask token>
print('length of fmri after removing skip:', len(initial_fmri_x))
print('fmri time runs from ', initial_fmri_x[0], ' to ', initial_fmri_x[-1])
print('reading in source files')
<mask token>
if not tide_io.checkspacedimmatch(immeandims, fmridims):
    print('immean file does not match')
    exit()
<mask token>
if not tide_io.checkspacedimmatch(boldpcdims, fmridims):
    print('boldpc file does not match')
    exit()
<mask token>
if not tide_io.checkspacedimmatch(lagdims, fmridims):
    print('lag file does not match')
    exit()
<mask token>
for k in range(0, numslices):
    sliceoffsettime = tide_resample.calcsliceoffset(sliceordertype, k,
        numslices, fmritr)
    print('processing slice ', k, ': sliceoffsettime=', sliceoffsettime)
    for j in range(0, ysize):
        for i in range(0, xsize):
            fmri_x = initial_fmri_x - lagdata[i, j, k] - sliceoffsettime
            if debug:
                print(i, j, k, fmri_x[0], initial_fmri_x[0], lagdata[i, j,
                    k], sliceoffsettime)
            if fastresample:
                fmri_y = genlagtc.yfromx(fmri_x)
            else:
                fmri_y = tide_resample.doresample(nirs_x, nirs_y, fmri_x)
            thenoise = noiselevel * standard_normal(len(fmri_y))
            simdata[i, j, k, :] = immeandata[i, j, k] * (1.0 + boldpcdata[i,
                j, k] / 100.0 * fmri_y) + thenoise
tide_io.savetonifti(simdata, fmriheader, fmrisizes, outputname)

<mask token>


def usage():
    print(
        'usage: simdata fmrifilename immeanfilename boldpcfilename lagfilename regressor regressorsamprate numskip inputstarttime outputname noiselevel sliceorder'
        )
    print('')
    print('required arguments:')
    print('\tfmrifilename\t\t- the example BOLD fmri file')
    print('\timmeanfilename\t\t- the 3D mean image file')
    print('\tboldpcfilename\t\t- the 3D bold percentage file')
    print('\tlagfilename\t\t- the 3D time delay (in seconds) at every point')
    print('\tregressorfile\t\t- the name of the regressor text file')
    print('\tregressorsamprate\t- the sample rate of the regressor file, in Hz'
        )
    print(
        '\tnumskip\t\t\t- the use to simulate tr periods deleted during preprocessing'
        )
    print(
        '\tinputstarttime\t\t- the time delay, in seconds, into the inputfile that matches the start time of the fmrifile'
        )
    print('\toutputname\t\t- the root name for the output files')
    print('\tnoiselevel\t\t- the variance of the noise')
    print('\tsliceorder\t\t- slice acquisition order:')
    print('\t\t\t\t\t0 : None')
    print('\t\t\t\t\t1 : Regular up (0, 1, 2, 3, ...)')
    print('\t\t\t\t\t2 : Regular down')
    print('\t\t\t\t\t3 : Use slice order file')
    print('\t\t\t\t\t4 : Use slice timings file')
    print('\t\t\t\t\t5 : Standard Interleaved (0, 2, 4 ... 1, 3, 5 ... )')
    print(
        '\t\t\t\t\t6 : Siemens Interleaved (0, 2, 4 ... 1, 3, 5 ... for odd number of slices)'
        )
    print(
        '\t\t\t\t\t                        (1, 3, 5 ... 0, 2, 4 ... for even number of slices)'
        )
    return ()


displayplots = False
fastresample = True
debug = False
if len(sys.argv) != 12:
    usage()
    exit()
fmrifilename = sys.argv[1]
immeanfilename = sys.argv[2]
boldpcfilename = sys.argv[3]
lagfilename = sys.argv[4]
filename = sys.argv[5]
inputfreq = float(sys.argv[6])
fmriskip = int(sys.argv[7])
inputstarttime = float(sys.argv[8])
outputname = sys.argv[9]
noiselevel = float(sys.argv[10])
sliceordertype = int(sys.argv[11])
print('fmrifilename', fmrifilename)
print('immeanfilename', immeanfilename)
print('boldpcfilename', boldpcfilename)
print('lagfilename', lagfilename)
print('filename', filename)
print('inputfreq', inputfreq)
print('fmriskip', fmriskip)
print('inputstarttime', inputstarttime)
print('outputname', outputname)
print('noiselevel', noiselevel)
fmritr, numtrs = tide_io.fmritimeinfo(fmrifilename)
nim_fmri, fmridata, fmriheader, fmridims, fmrisizes = tide_io.readfromnifti(
    fmrifilename)
print('fmri data: ', numtrs, ' timepoints, tr = ', fmritr)
inputvec = tide_math.stdnormalize(tide_io.readvec(filename))
simregressorpts = len(inputvec)
print('Input regressor has ', simregressorpts, ' points')
inputstep = 1.0 / inputfreq
nirs_x = r_[0.0:1.0 * simregressorpts] * inputstep - inputstarttime
nirs_y = inputvec[0:simregressorpts]
print('nirs regressor runs from ', nirs_x[0], ' to ', nirs_x[-1])
fmrifreq = 1.0 / fmritr
initial_fmri_x = r_[0.0:fmritr * (numtrs - fmriskip):fmritr
    ] + fmritr * fmriskip
print('length of fmri after removing skip:', len(initial_fmri_x))
print('fmri time runs from ', initial_fmri_x[0], ' to ', initial_fmri_x[-1])
print('reading in source files')
nim_immean, immeandata, immeanheader, immeandims, immeansizes = (tide_io.
    readfromnifti(immeanfilename))
if not tide_io.checkspacedimmatch(immeandims, fmridims):
    print('immean file does not match')
    exit()
nim_boldpc, boldpcdata, boldpcheader, boldpcdims, boldpcsizes = (tide_io.
    readfromnifti(boldpcfilename))
if not tide_io.checkspacedimmatch(boldpcdims, fmridims):
    print('boldpc file does not match')
    exit()
nim_lag, lagdata, lagheader, lagdims, lagsizes = tide_io.readfromnifti(
    lagfilename)
if not tide_io.checkspacedimmatch(lagdims, fmridims):
    print('lag file does not match')
    exit()
thesizes = fmrisizes
xdim = thesizes[0]
ydim = thesizes[1]
slicethickness = thesizes[2]
tr = thesizes[3]
thedims = fmridims
xsize = thedims[1]
ysize = thedims[2]
numslices = thedims[3]
timepoints = thedims[4]
simdata = zeros((xsize, ysize, numslices, len(initial_fmri_x)), dtype='float')
fmrilcut = 0.0
fmriucut = fmrifreq / 2.0
padvalue = 60.0
numpadtrs = int(padvalue / fmritr)
padvalue = fmritr * numpadtrs
genlagtc = tide_resample.fastresampler(nirs_x, nirs_y, padvalue=padvalue,
    doplot=False)
fmri_y_f = genlagtc.yfromx(initial_fmri_x)
fmri_y = tide_resample.doresample(nirs_x, nirs_y, initial_fmri_x)
for k in range(0, numslices):
    sliceoffsettime = tide_resample.calcsliceoffset(sliceordertype, k,
        numslices, fmritr)
    print('processing slice ', k, ': sliceoffsettime=', sliceoffsettime)
    for j in range(0, ysize):
        for i in range(0, xsize):
            fmri_x = initial_fmri_x - lagdata[i, j, k] - sliceoffsettime
            if debug:
                print(i, j, k, fmri_x[0], initial_fmri_x[0], lagdata[i, j,
                    k], sliceoffsettime)
            if fastresample:
                fmri_y = genlagtc.yfromx(fmri_x)
            else:
                fmri_y = tide_resample.doresample(nirs_x, nirs_y, fmri_x)
            thenoise = noiselevel * standard_normal(len(fmri_y))
            simdata[i, j, k, :] = immeandata[i, j, k] * (1.0 + boldpcdata[i,
                j, k] / 100.0 * fmri_y) + thenoise
tide_io.savetonifti(simdata, fmriheader, fmrisizes, outputname)

from __future__ import print_function
import sys
import rapidtide.miscmath as tide_math
import rapidtide.io as tide_io
import rapidtide.resample as tide_resample
from matplotlib.pyplot import *


def usage():
    print(
        'usage: simdata fmrifilename immeanfilename boldpcfilename lagfilename regressor regressorsamprate numskip inputstarttime outputname noiselevel sliceorder'
        )
    print('')
    print('required arguments:')
    print('\tfmrifilename\t\t- the example BOLD fmri file')
    print('\timmeanfilename\t\t- the 3D mean image file')
    print('\tboldpcfilename\t\t- the 3D bold percentage file')
    print('\tlagfilename\t\t- the 3D time delay (in seconds) at every point')
    print('\tregressorfile\t\t- the name of the regressor text file')
    print('\tregressorsamprate\t- the sample rate of the regressor file, in Hz'
        )
    print(
        '\tnumskip\t\t\t- the use to simulate tr periods deleted during preprocessing'
        )
    print(
        '\tinputstarttime\t\t- the time delay, in seconds, into the inputfile that matches the start time of the fmrifile'
        )
    print('\toutputname\t\t- the root name for the output files')
    print('\tnoiselevel\t\t- the variance of the noise')
    print('\tsliceorder\t\t- slice acquisition order:')
    print('\t\t\t\t\t0 : None')
    print('\t\t\t\t\t1 : Regular up (0, 1, 2, 3, ...)')
    print('\t\t\t\t\t2 : Regular down')
    print('\t\t\t\t\t3 : Use slice order file')
    print('\t\t\t\t\t4 : Use slice timings file')
    print('\t\t\t\t\t5 : Standard Interleaved (0, 2, 4 ... 1, 3, 5 ... )')
    print(
        '\t\t\t\t\t6 : Siemens Interleaved (0, 2, 4 ... 1, 3, 5 ... for odd number of slices)'
        )
    print(
        '\t\t\t\t\t                        (1, 3, 5 ... 0, 2, 4 ... for even number of slices)'
        )
    return ()


displayplots = False
fastresample = True
debug = False
if len(sys.argv) != 12:
    usage()
    exit()
fmrifilename = sys.argv[1]
immeanfilename = sys.argv[2]
boldpcfilename = sys.argv[3]
lagfilename = sys.argv[4]
filename = sys.argv[5]
inputfreq = float(sys.argv[6])
fmriskip = int(sys.argv[7])
inputstarttime = float(sys.argv[8])
outputname = sys.argv[9]
noiselevel = float(sys.argv[10])
sliceordertype = int(sys.argv[11])
print('fmrifilename', fmrifilename)
print('immeanfilename', immeanfilename)
print('boldpcfilename', boldpcfilename)
print('lagfilename', lagfilename)
print('filename', filename)
print('inputfreq', inputfreq)
print('fmriskip', fmriskip)
print('inputstarttime', inputstarttime)
print('outputname', outputname)
print('noiselevel', noiselevel)
fmritr, numtrs = tide_io.fmritimeinfo(fmrifilename)
nim_fmri, fmridata, fmriheader, fmridims, fmrisizes = tide_io.readfromnifti(
    fmrifilename)
print('fmri data: ', numtrs, ' timepoints, tr = ', fmritr)
inputvec = tide_math.stdnormalize(tide_io.readvec(filename))
simregressorpts = len(inputvec)
print('Input regressor has ', simregressorpts, ' points')
inputstep = 1.0 / inputfreq
nirs_x = r_[0.0:1.0 * simregressorpts] * inputstep - inputstarttime
nirs_y = inputvec[0:simregressorpts]
print('nirs regressor runs from ', nirs_x[0], ' to ', nirs_x[-1])
fmrifreq = 1.0 / fmritr
initial_fmri_x = r_[0.0:fmritr * (numtrs - fmriskip):fmritr
    ] + fmritr * fmriskip
print('length of fmri after removing skip:', len(initial_fmri_x))
print('fmri time runs from ', initial_fmri_x[0], ' to ', initial_fmri_x[-1])
print('reading in source files')
nim_immean, immeandata, immeanheader, immeandims, immeansizes = (tide_io.
    readfromnifti(immeanfilename))
if not tide_io.checkspacedimmatch(immeandims, fmridims):
    print('immean file does not match')
    exit()
nim_boldpc, boldpcdata, boldpcheader, boldpcdims, boldpcsizes = (tide_io.
    readfromnifti(boldpcfilename))
if not tide_io.checkspacedimmatch(boldpcdims, fmridims):
    print('boldpc file does not match')
    exit()
nim_lag, lagdata, lagheader, lagdims, lagsizes = tide_io.readfromnifti(
    lagfilename)
if not tide_io.checkspacedimmatch(lagdims, fmridims):
    print('lag file does not match')
    exit()
thesizes = fmrisizes
xdim = thesizes[0]
ydim = thesizes[1]
slicethickness = thesizes[2]
tr = thesizes[3]
thedims = fmridims
xsize = thedims[1]
ysize = thedims[2]
numslices = thedims[3]
timepoints = thedims[4]
simdata = zeros((xsize, ysize, numslices, len(initial_fmri_x)), dtype='float')
fmrilcut = 0.0
fmriucut = fmrifreq / 2.0
padvalue = 60.0
numpadtrs = int(padvalue / fmritr)
padvalue = fmritr * numpadtrs
genlagtc = tide_resample.fastresampler(nirs_x, nirs_y, padvalue=padvalue,
    doplot=False)
fmri_y_f = genlagtc.yfromx(initial_fmri_x)
fmri_y = tide_resample.doresample(nirs_x, nirs_y, initial_fmri_x)
for k in range(0, numslices):
    sliceoffsettime = tide_resample.calcsliceoffset(sliceordertype, k,
        numslices, fmritr)
    print('processing slice ', k, ': sliceoffsettime=', sliceoffsettime)
    for j in range(0, ysize):
        for i in range(0, xsize):
            fmri_x = initial_fmri_x - lagdata[i, j, k] - sliceoffsettime
            if debug:
                print(i, j, k, fmri_x[0], initial_fmri_x[0], lagdata[i, j,
                    k], sliceoffsettime)
            if fastresample:
                fmri_y = genlagtc.yfromx(fmri_x)
            else:
                fmri_y = tide_resample.doresample(nirs_x, nirs_y, fmri_x)
            thenoise = noiselevel * standard_normal(len(fmri_y))
            simdata[i, j, k, :] = immeandata[i, j, k] * (1.0 + boldpcdata[i,
                j, k] / 100.0 * fmri_y) + thenoise
tide_io.savetonifti(simdata, fmriheader, fmrisizes, outputname)

#!/usr/bin/env python
#
#   Copyright 2016 Blaise Frederick
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.
#
#
#       $Author: frederic $
#       $Date: 2016/06/14 12:04:51 $
#       $Id: simdata,v 1.18 2016/06/14 12:04:51 frederic Exp $
#
from __future__ import print_function
import sys
import rapidtide.miscmath as tide_math
import rapidtide.io as tide_io
import rapidtide.resample as tide_resample
from matplotlib.pyplot import *


def usage():
    print(
        "usage: simdata fmrifilename immeanfilename boldpcfilename lagfilename regressor regressorsamprate numskip inputstarttime outputname noiselevel sliceorder")
    print("")
    print("required arguments:")
    print("	fmrifilename		- the example BOLD fmri file")
    print("	immeanfilename		- the 3D mean image file")
    print("	boldpcfilename		- the 3D bold percentage file")
    print("	lagfilename		- the 3D time delay (in seconds) at every point")
    print("	regressorfile		- the name of the regressor text file")
    print("	regressorsamprate	- the sample rate of the regressor file, in Hz")
    print("	numskip			- the use to simulate tr periods deleted during preprocessing")
    print(
        "	inputstarttime		- the time delay, in seconds, into the inputfile that matches the start time of the fmrifile")
    print("	outputname		- the root name for the output files")
    print("	noiselevel		- the variance of the noise")
    print("	sliceorder		- slice acquisition order:")
    print("					0 : None")
    print("					1 : Regular up (0, 1, 2, 3, ...)")
    print("					2 : Regular down")
    print("					3 : Use slice order file")
    print("					4 : Use slice timings file")
    print("					5 : Standard Interleaved (0, 2, 4 ... 1, 3, 5 ... )")
    print("					6 : Siemens Interleaved (0, 2, 4 ... 1, 3, 5 ... for odd number of slices)")
    print("					                        (1, 3, 5 ... 0, 2, 4 ... for even number of slices)")

    return ()


# set default variable values
displayplots = False
fastresample = True
debug = False

# get the command line parameters
if len(sys.argv) != 12:
    usage()
    exit()

# handle required args first
fmrifilename = sys.argv[1]
immeanfilename = sys.argv[2]
boldpcfilename = sys.argv[3]
lagfilename = sys.argv[4]
filename = sys.argv[5]
inputfreq = float(sys.argv[6])
fmriskip = int(sys.argv[7])
inputstarttime = float(sys.argv[8])
outputname = sys.argv[9]
noiselevel = float(sys.argv[10])
sliceordertype = int(sys.argv[11])

print("fmrifilename", fmrifilename)
print("immeanfilename", immeanfilename)
print("boldpcfilename", boldpcfilename)
print("lagfilename", lagfilename)
print("filename", filename)
print("inputfreq", inputfreq)
print("fmriskip", fmriskip)
print("inputstarttime", inputstarttime)
print("outputname", outputname)
print("noiselevel", noiselevel)

fmritr, numtrs = tide_io.fmritimeinfo(fmrifilename)
nim_fmri, fmridata, fmriheader, fmridims, fmrisizes = tide_io.readfromnifti(fmrifilename)
print("fmri data: ", numtrs, " timepoints, tr = ", fmritr)

# read in the timecourse to resample
inputvec = tide_math.stdnormalize(tide_io.readvec(filename))
simregressorpts = len(inputvec)

# prepare the input data for interpolation
print("Input regressor has ", simregressorpts, " points")
inputstep = 1.0 / inputfreq
nirs_x = r_[0.0:1.0 * simregressorpts] * inputstep - inputstarttime
nirs_y = inputvec[0:simregressorpts]
print('nirs regressor runs from ', nirs_x[0], ' to ', nirs_x[-1])

# prepare the output timepoints
fmrifreq = 1.0 / fmritr
initial_fmri_x = r_[0.0:fmritr * (numtrs - fmriskip):fmritr] + fmritr * fmriskip
print('length of fmri after removing skip:', len(initial_fmri_x))
print('fmri time runs from ', initial_fmri_x[0], ' to ', initial_fmri_x[-1])

# read in the immean, boldpc and lag files
print("reading in source files")
nim_immean, immeandata, immeanheader, immeandims, immeansizes = tide_io.readfromnifti(immeanfilename)
if not tide_io.checkspacedimmatch(immeandims, fmridims):
    print("immean file does not match")
    exit()
nim_boldpc, boldpcdata, boldpcheader, boldpcdims, boldpcsizes = tide_io.readfromnifti(boldpcfilename)
if not tide_io.checkspacedimmatch(boldpcdims, fmridims):
    print("boldpc file does not match")
    exit()
nim_lag, lagdata, lagheader, lagdims, lagsizes = tide_io.readfromnifti(lagfilename)
if not tide_io.checkspacedimmatch(lagdims, fmridims):
    print("lag file does not match")
    exit()

# now set up the simulated data array
thesizes = fmrisizes
xdim = thesizes[0]
ydim = thesizes[1]
slicethickness = thesizes[2]
tr = thesizes[3]

thedims = fmridims
xsize = thedims[1]
ysize = thedims[2]
numslices = thedims[3]
timepoints = thedims[4]

simdata = zeros((xsize, ysize, numslices, len(initial_fmri_x)), dtype='float')

fmrilcut = 0.0
fmriucut = fmrifreq / 2.0

# set up fast resampling
padvalue = 60.0
numpadtrs = int(padvalue / fmritr)
padvalue = fmritr * numpadtrs

genlagtc = tide_resample.fastresampler(nirs_x, nirs_y, padvalue=padvalue, doplot=False)
fmri_y_f = genlagtc.yfromx(initial_fmri_x)
fmri_y = tide_resample.doresample(nirs_x, nirs_y, initial_fmri_x)
# fig = figure()
# ax = fig.add_subplot(111)
# plot(initial_fmri_x, fmri_y_f + 1.0, initial_fmri_x, fmri_y)
# show()

# fig = figure()
# ax = fig.add_subplot(111)
# ax.set_title('Regressors')
# plot(nirs_x, nirs_y)
# show()

# loop over space
for k in range(0, numslices):
    sliceoffsettime = tide_resample.calcsliceoffset(sliceordertype, k, numslices, fmritr)
    print("processing slice ", k, ": sliceoffsettime=", sliceoffsettime)
    for j in range(0, ysize):
        for i in range(0, xsize):
            # correct the output time points
            fmri_x = initial_fmri_x - lagdata[i, j, k] - sliceoffsettime
            if debug:
                print(i, j, k, fmri_x[0], initial_fmri_x[0], lagdata[i, j, k], sliceoffsettime)

            if fastresample:
                fmri_y = genlagtc.yfromx(fmri_x)
            else:
                fmri_y = tide_resample.doresample(nirs_x, nirs_y, fmri_x)
            thenoise = noiselevel * standard_normal(len(fmri_y))
            simdata[i, j, k, :] = immeandata[i, j, k] * (1.0 + (boldpcdata[i, j, k] / 100.0) * fmri_y) + thenoise

tide_io.savetonifti(simdata, fmriheader, fmrisizes, outputname)
