Ls user function

In [ ]:

import math
import numpy as np
import pint
ureg = pint.UnitRegistry()
from dynasaur.calc.cfc import CFC
from scipy.integrate import cumtrapz
from scipy.stats import norm 
from scipy import special

import math import numpy as np import pint ureg = pint.UnitRegistry() from dynasaur.calc.cfc import CFC from scipy.integrate import cumtrapz from scipy.stats import norm from scipy import special

In [ ]:

from dynasaur.calc.object_calculation_util import ObjectCalcUtil, UniversalLimit
from dynasaur.utils.constants import LOGConstants, PluginsParamDictDef, OutputStringForPlugins

from dynasaur.calc.object_calculation_util import ObjectCalcUtil, UniversalLimit from dynasaur.utils.constants import LOGConstants, PluginsParamDictDef, OutputStringForPlugins

In [ ]:

# use to define own functions
class UserFunction:



    @staticmethod
    def time_at_Xmove(param_dict, units):
        move = (param_dict['movement'])
        time = (param_dict['time'])
        value = (param_dict['value'])
        offset = (param_dict['offset'])

        timestep_of_maximum = 0
        for t in range(len(move)):
            y = move[t]
            if y > value:
               timestep_of_maximum = t
               break
               
        time_of_maximum = time[timestep_of_maximum] - offset
       
        return (np.max(time_of_maximum))

    @staticmethod
    def vel_at_timeX(param_dict, units):
        move = (param_dict['movement'])
        vel = (param_dict['velocity'])
        value = (param_dict['value'])

        for t in range(len(move)):
            y = move[t]
            if y > value:
                timestep_of_maximum = t
                break

        velocity_at_timeX = vel[timestep_of_maximum]

        return (np.max(velocity_at_timeX))

    @staticmethod
    def FC_t(param_dict, units):

        t = param_dict['time'].flatten()  # / units.second()  # t in [s]
        fr = param_dict['Fres'].flatten()
        fbmax = max(fr)

        f_rise = fbmax * 0.01
        trise = -1

        for i in range(len(t) - 1):
            if f_rise < fr[i]:
                trise = t[i]
                break

        # inclusive 20ms start_time! -- zeitpunkt in der fbelt um 20% gefallen ist
        return trise


    @staticmethod
    def H350_FEMUR_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        # FEMUR scalar value - injuryvalue
        risk = 1 / (1 + math.exp(offset - scaling * (np.abs(injuryvalue))))

        return (risk)


    @staticmethod
    def H350_NECKCOMP_AIS3(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * (np.abs(injuryvalue[0]))))

        return (risk)


    @staticmethod
    def H350_NECKTEN_AIS3(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * (np.abs(injuryvalue[0]))))

        return (risk)


    @staticmethod
    def H350_CHESTDEF_AIS3(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * ((np.abs(injuryvalue[0])) ** power)))

        return (risk)


    @staticmethod
    def H350_NIJ_AIS3(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]

        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * injuryvalue[0]))

        return (risk)


    @staticmethod
    def THOR50_NIJ_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * injuryvalue[0]))

        return (risk)


    @staticmethod
    def THOR50_NIJ_AIS3(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]

        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * injuryvalue[0]))

        # "name" : "THOR50_NIJ_AIS3+",
        # "ccdf" : "1 / (1 + math.exp(4.9372 - 4.5294*value))"

        return (risk)


    @staticmethod
    def THOR50_CHESTDEF_AIS3(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 - math.exp(-((injuryvalue / divide) ** power))

        # "name" : "THOR50_CHESTDEF_AIS3+",
        # "ccdf" : "1 - math.exp(-((value/59.865)**2.7187))"

        return (risk)


    @staticmethod
    def THOR50_ABDOMENDEF_AIS3(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]

        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * injuryvalue))

        # "name" : "THOR50_ABDOMENDEF_AIS3+",
        # "ccdf" : "1 / (1 + math.exp(7.849 - 0.0886 * value))"

        return (risk)


    @staticmethod
    def THOR50_FEMUR_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * (abs(injuryvalue))))

        # "name" : "THOR50_FEMUR_AIS2+",
        # "ccdf" : "1 / (1 + math.exp(5.7949 - 0.6748*(abs(value))))"

        return (risk)


    @staticmethod
    def THOR50_TIBIA_FZ_UPPER_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]

        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * injuryvalue))

        # "name" : "THOR50_TIBIA_FZ_UPPER_AIS2+",
        # "ccdf" : "1 / (1 + math.exp(5.6654 - 0.8189*value))"

        return (risk)


    @staticmethod
    def THOR50_TIBIA_FZ_LOWER_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 / (1 + math.exp(offset - scaling * injuryvalue))

        # "name" : "THOR50_TIBIA_FZ_LOWER_AIS2+",
        # "ccdf" : "1 / (1 + math.exp(3.9121 - 0.48*value))"

        return (risk)


    @staticmethod
    def THOR50_TIBIA_Mres_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = 1 - math.exp((math.log(injuryvalue) - offset) / divide)

        # "name" : "THOR50_TIBIA_Mres_AIS2+",
        # "ccdf" : "1 - math.exp( - math.exp((math.log(value)-5.8492)/0.2965))"

        return (risk)


    @staticmethod
    def THOR50_HIP_Fracture_AIS(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]

        risk = 0
        risk = (1.0 + math.erf(((math.log(scaling * injuryvalue[0]) - offset) / divide) / math.sqrt(2.0))) / 2.0

        # "name" : "THOR50_HIP_Fracture_AIS+",
        # "ccdf" : "(1.0 + math.erf(((math.log(1.429*value) - 1.6058) / 0.2339) / math.sqrt(2.0))) / 2.0"

        return (risk)


    @staticmethod
    def WS50_SHOULDERFORCE_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]
        age = param_dict["age"]

        risk = 0
        risk = 1 - math.exp(-(np.abs(injuryvalue) / (offset - scaling * age)) ** power)

        # "name" : "WS50_SHOULDERFORCE_AIS2+",
        # "ccdf" : "1 - (math.exp(-(value / (8.144 - 0.006 * 50)) ** 7.41))"

        return (risk)


    @staticmethod
    def WS50_ABDOMENDEF_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]
        age = param_dict["age"]

        risk = 0
        risk = 1 - math.exp(-((np.abs(injuryvalue[0]) / (offset - scaling * age))) ** power)
        # "name" : "WS50_ABDOMENDEF_AIS2+" Soft Tissue Abdominal Injury, 
        # "ccdf" : "1 - math.exp( - (injuryvalue[0] / (5.368 - 0.021 * 50)) ** 8.61)"
        return (risk)


    @staticmethod
    def WS50_DEF_SKEL_AIS3(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]
        age = param_dict["age"]

        risk = 0
        risk = 1 / (1 + math.exp( -(math.log(abs(injuryvalue[0])) - (offset - scaling * age)) / divide))
        
        
        return (risk)

    @staticmethod
    def WS50_PUBICFORCE_AIS2(param_dict, units):
        injuryvalue = param_dict["injuryvalue"]
        offset = param_dict["offset"]
        divide = param_dict["divide"]
        scaling = param_dict["scaling"]
        power = param_dict["power"]
        age = param_dict["age"]

        risk = 0
        risk = 1 - math.exp( - ( ( injuryvalue[0] / (offset - scaling * age) ) ** power) )

        # "name" : "WS50_PUBICFORCE_AIS2+",
        # "ccdf" : "1 - math.exp(-(value / (8.775 - 0.014 * 50)) ** 4.60)"

        return (risk)

    @staticmethod
    def THOR_TCI_Rmax(Rmax_Points_UL: pint.Quantity, Rmax_Points_LL: pint.Quantity, Rmax_Points_UR: pint.Quantity, Rmax_Points_LR: pint.Quantity):

        Rmax = 0 * Rmax_Points_UL.units
        
        Rmax = max([np.max(Rmax_Points_UL[:, 0]),
                    np.max(Rmax_Points_LL[:, 0]),
                    np.max(Rmax_Points_UR[:, 0]),
                    np.max(Rmax_Points_LR[:, 0])])
        return (Rmax)

    @staticmethod
    def THOR_Abdomen_Deflection(Abdmax_LEFT: pint.Quantity, Abdmax_RIGHT: pint.Quantity):

        Abdomenmax = 0 * Abdmax_LEFT.units
        
        Abdomenmax = max([np.max(Abdmax_RIGHT[:, 0]),
                          np.max(Abdmax_LEFT[:, 0])])

        return (Abdomenmax)

    @staticmethod
    def AIRBAG_CONTACT(param_dict, units):
        data_vector = param_dict['data_vector']
        time = param_dict['time']
        move = param_dict['move']
        tbase = param_dict['tbase']

        t_cont = 0
        for i in range(len(time) - 1):
            if move < np.abs(data_vector[i]):
                t_cont = time[i] - tbase
                break

        return np.max(t_cont)

    @staticmethod
    def VEL_NODE_MAX(param_dict, units):

        t = param_dict['time'].flatten()
        y = param_dict['y'].flatten()

        vel = []
        velmax = 0
        vel.append(velmax)

        for i in range(len(t) - 1):
            if i > 0:
                vel.append((y[i] - y[i - 1]) / (t[i] - t[i - 1]))
            if np.abs(vel[i]) > np.abs(velmax):
                velmax = vel[i]

        return (velmax)

    @staticmethod
    def VEL_NODE(param_dict, units):

        t = param_dict['time'].flatten()
        y = param_dict['y'].flatten()

        vel = []
        velmax = 0
        vel.append(velmax)

        for i in range(len(t) - 1):
            if i > 0:
                vel.append((y[i] - y[i - 1]) / (t[i] - t[i - 1]))

        return np.array(vel).reshape(-1, 1)

    @staticmethod
    def ext_bend_mom(param_dict, units):
        force_x = param_dict['force_x'] / (units.meter() / (units.second() ** 2))
        moment_y = param_dict['moment_y'] / ((units.meter() ** 2) / (units.second() ** 2))
        moc_d = param_dict['distance_occipital_condyle'] / ((units.meter() ** 2) / (units.second() ** 2))
        mtot = moment_y - moc_d * force_x
        return (mtot)

    @staticmethod
    def double_abs_sub_res(arr1: pint.Quantity, arr2: pint.Quantity, arr3: pint.Quantity, arr4: float):          
        arr1m = arr1.magnitude
        arr2m = arr2.magnitude
        arr3m = arr3.magnitude
        arr4m = arr4.magnitude
        length = arr3.units
        
        rel_dist_max_1 = []
        rel_dist_max_2 = []
        rel_dist = []

        rel_dist_1 = np.subtract(arr1m, arr2m)
        rel_dist_2 = np.subtract(arr1m, arr3m)
                
        for i in range(len(rel_dist_1)):
            rel_dist_max_1.append(
                (rel_dist_1[i, 0] ** 2 + rel_dist_1[i, 1] ** 2 + rel_dist_1[i, 2] ** 2) ** (1 / 2) - arr4m)
        for i in range(len(rel_dist_2)):
            rel_dist_max_2.append(
                (rel_dist_2[i, 0] ** 2 + rel_dist_2[i, 1] ** 2 + rel_dist_2[i, 2] ** 2) ** (1 / 2) - arr4m)
        if np.min(rel_dist_max_1) < 0 or np.min(rel_dist_max_2) < 0:
            for i in range(len(rel_dist_1)):
                rel_dist.append(0)
                                        
        else:
            if np.max(np.abs(rel_dist_max_1)) < np.max(np.abs(rel_dist_max_2)):
                rel_dist = rel_dist_max_1
            else:
                rel_dist = rel_dist_max_2                 
        
        
        return np.abs(rel_dist) * length

    @staticmethod
    def MOCy_max(param_dict, units):
        force_x = param_dict['force_x'] / (units.meter() / (units.second() ** 2))
        moment_y = param_dict['moment_y'] / ((units.meter() ** 2) / (units.second() ** 2))
        moc_d = param_dict['distance_occipital_condyle']
        MOCy_max = np.max(moment_y - moc_d * force_x)
        return (MOCy_max)

    @staticmethod
    def MOCy_min(param_dict, units):
        force_x = param_dict['force_x'] / (units.meter() / (units.second() ** 2))
        moment_y = param_dict['moment_y'] / ((units.meter() ** 2) / (units.second() ** 2))
        moc_d = param_dict['distance_occipital_condyle']
        MOCy_min = np.min(moment_y - moc_d * force_x)
        return (MOCy_min)

    @staticmethod
    def minimum_of_four(data_vector_1: pint.Quantity, data_vector_2: pint.Quantity, data_vector_3: pint.Quantity, data_vector_4: pint.Quantity):
        min_1 = (np.min(data_vector_1)).flatten()
        min_2 = (np.min(data_vector_2)).flatten()
        min_3 = (np.min(data_vector_3)).flatten()
        min_4 = (np.min(data_vector_4)).flatten()
        temp = [min_1, min_2, min_3, min_4]
        min = np.min(temp)
        # index_of_minimum = np.argmin(param_dict['minimum'])
        # time_of_minimum = param_dict['time'][index_of_minimum]

        return (min) * data_vector_1.units

    @staticmethod
    def maximum_of_four(data_vector_1: pint.Quantity, data_vector_2: pint.Quantity, data_vector_3: pint.Quantity, data_vector_4: pint.Quantity):
        max_1 = (np.max(data_vector_1)).flatten()
        max_2 = (np.max(data_vector_2)).flatten()
        max_3 = (np.max(data_vector_3)).flatten()
        max_4 = (np.max(data_vector_4)).flatten()
        temp = [max_1, max_2, max_3, max_4]
        max = np.max(temp)
        # index_of_minimum = np.argmin(param_dict['minimum'])
        # time_of_minimum = param_dict['time'][index_of_minimum]

        return (max) * data_vector_1.units

    @staticmethod
    def ABS_MAX_VALUE(param_dict, units):
        data_vec = param_dict['data_vector']
        val = 0
        if np.abs(np.min(data_vec)) < np.abs(np.max(data_vec)):
            val = np.max(data_vec)
        else:
            val = np.min(data_vec)
        return val

    @staticmethod
    def RES_ACC_from_XYZ_Velocities(a_x: pint.Quantity, a_y: pint.Quantity, a_z: pint.Quantity):    
        a_x = param_dict['a_x'].flatten()
        a_y = param_dict['a_y'].flatten()
        a_z = param_dict['a_z'].flatten()
        a_res = (a_x ** 2 + a_y ** 2 + a_z ** 2) ** (1 / 2)

        return np.array(a_res).reshape(-1, 1)

    @staticmethod
    def ILIAC_Drop(Fz: pint.Quantity, time: pint.Quantity):
           
        t = time.to_base_units().magnitude 
        fi = Fz.magnitude    
        iliac_drop = []
        delta_t = (t[1] - t[0])[0]        
        
        tresample = int(round(0.001 / delta_t))
        for i in range(len(t) - (1 + tresample)):
            value = fi[i + tresample] - fi[i]
            if delta_t < 0.001 and value < 0:
                iliac_drop.append(value[0])
            else:
                iliac_drop.append(0)
        
        return np.array(iliac_drop).reshape(-1, 1) * time.units



    @staticmethod
    def FB_drop_t(time: pint.Quantity, Fb: pint.Quantity, drop_value_percent: pint.Quantity):
    
        t = time
        fb = Fb
        dvp = drop_value_percent

        fbmax = max(fb)
        fbmax_drop = fbmax * (1 - dvp / 100)
        tswitch = 0
        fb_avg = 0

        for i in range(len(t) - 1):

            fb_avg = fb_avg + fb[i]
            if fbmax == fb[i]:
                tmax = t[i]
                tswitch = 1

            if tswitch == 1:
                if fbmax_drop >= fb[i]:
                    tdrop = t[i]
                    average_i = i
                    tswitch = 0
                    break
            if tswitch == 0:
                tdrop = t[i]
                average_i = i

        # durchschnittliche Kraft bis fbmax_drop
        fb_avg_drop = fb_avg / average_i

        # inclusive 20ms start_time! -- zeitpunkt in der fbelt um 20% gefallen ist
        return tdrop

    @staticmethod
    def FB_drop_avg(Fb: pint.Quantity, time: pint.Quantity, drop_value_percent: pint.Quantity):
 
        t = time
        fb = Fb
        dvp = drop_value_percent.magnitude

        fbmax = max(fb)

        fbmax_drop = fbmax * (1 - dvp / 100)
        tswitch = 0
        fb_avg = 0

        for i in range(len(t) - 1):

            fb_avg = fb_avg + fb[i]
            if fbmax == fb[i]:
                tmax = t[i]
                tswitch = 1

            if tswitch == 1:
                if fbmax_drop >= fb[i]:
                    tdrop = t[i]
                    average_i = i
                    tswitch = 0
                    break

            if tswitch == 0:
                tdrop = t[i]
                average_i = i

        # durchschnittliche Kraft bis fbmax_drop
        fb_avg_drop = fb_avg / average_i
        # nur ein return
        return fb_avg_drop

        # END Peter L. - BeltLoad80%###########
                
    @staticmethod
    def HIC_AIS2(HIC: pint.Quantity):
        HICAIS2 = norm.cdf((math.log(HIC) - 6.96352) / 0.84664)
        return HICAIS2       

    @staticmethod
    def HIC_AIS3(HIC: pint.Quantity):
        HICAIS3 = norm.cdf((math.log(HIC) - 7.45231) / 0.73998)
        return HICAIS3             
  
    @staticmethod
    def HIC_AIS4(HIC: pint.Quantity):
        HICAIS4 = norm.cdf((math.log(HIC) - 7.65605) / 0.60580)
        return HICAIS4   

    @staticmethod
    def BRIC_AIS4(BRIC: pint.Quantity):
        BRICAIS4 = 1 - np.exp((BRIC-0.523)/0.647) ** 1.8
        return BRICAIS4   
        
    @staticmethod
    def NIJ_AIS2_THOR(NIJ: pint.Quantity):
        NIJAIS2 = 1 / (1 + np.exp(4.3085 - 5.4079 * NIJ))
        return NIJAIS2                        

    @staticmethod
    def NIJ_AIS3_THOR(NIJ: pint.Quantity):
        NIJAIS3 = 1 / (1 + np.exp(4.9372 - 4.5294 * NIJ))
        return NIJAIS3

    @staticmethod
    def NIJ_AIS3_H305(NIJ: pint.Quantity):
        NIJAIS3 = 1 - np.exp(-(NIJ/1.3933) ** 2.8816)
        return NIJAIS3
    
    @staticmethod
    def THOR_TCI_Rmax_AIS3(TCI_Rmax: pint.Quantity):
        TCI_RmaxAIS3 = 1 - np.exp(-(THOR_TCI_Rmax_AIS3/59.865) ** 2.7187)
        return TCI_RmaxAIS3        

    @staticmethod
    def Abdomen_Deflection_AIS3_THOR(Abdomen_Deflection: pint.Quantity):
        Abdomen_DeflectionAIS3 = 1 / (1 + np.exp(7.849 - 0.0886 * Abdomen_Deflection))
        return Abdomen_DeflectionAIS3   
        
    @staticmethod
    def Femur_axial_force_AIS2_THOR(Femur_Comp: pint.Quantity):
        Femur_Comp_mag = Femur_Comp.magnitude
        Femur_axialAIS2 = 1 / (1 + np.exp(5.7949 - 0.6748 * Femur_Comp_mag))
        return Femur_axialAIS2    
        
    @staticmethod
    def Femur_axial_force_AIS2_H305(Femur_Comp: pint.Quantity):
        Femur_Comp_mag = Femur_Comp.magnitude    
        Femur_axialAIS2 = 1 / (1 + np.exp(5.7949 - 0.7619 * Femur_Comp_mag))
        return Femur_axialAIS2          
                    
    @staticmethod
    def Lower_Tibia_axial_force_AIS2_THOR(Tibia_Fz: pint.Quantity):
        Tibia_Fz_mag = Tibia_Fz.magnitude    
        Tibia_FzAIS2 = 1 / (1 + np.exp(3.9121 - 0.48 * Tibia_Fz_mag))
        return Tibia_FzAIS2         

    @staticmethod
    def Upper_Tibia_axial_force_AIS2_THOR(Tibia_Fz: pint.Quantity):
        Tibia_Fz_mag = Tibia_Fz.magnitude        
        Tibia_FzAIS2 = 1 / (1 + np.exp(5.6654 - 0.8189 * Tibia_Fz_mag))
        return Tibia_FzAIS2  
        
    @staticmethod
    def Tibia_moment_AIS2_THOR(Tibia_Moment: pint.Quantity):
        Tibia_Momentm =Tibia_Moment.magnitude
        Tibia_MomentAIS2 = 1 - np.exp(-np.exp((math.log(Tibia_Momentm) - 5.8492)/0.2965))
        return Tibia_MomentAIS2         
        
    @staticmethod
    def initValue (arr1: pint.Quantity):   
        initialValue = arr1[0]              
        return initialValue
        
    @staticmethod
    def res_wo_unit(arr1: pint.Quantity, arr2: pint.Quantity, arr3: pint.Quantity):
        arr1m = arr1.magnitude
        arr2m = arr2.magnitude
        arr3m = arr3.magnitude
        unit = arr1.units

        return np.sqrt(arr1m**2 + arr2m**2 + arr3m**2) * unit        
        
    @staticmethod
    def percentilemin(object_data: dict, selection_tension_compression: str,
                   integration_point: str, percentile: float):
        """only one value possible to return

        Args:
          object_data: param selection_tension_compression:
          integration_point: param percentile:
          object_data: dict: 
          selection_tension_compression: str: 
          integration_point: str: 
          percentile: str: 
          units: 

        Returns:

        """

        part_ids = object_data["part_ids"]
        part_data = object_data["part_data"]
        part_idx = object_data["part_idx"]

        if len(part_ids) == 0:
            return np.NaN

        obj_calc = ObjectCalcUtil()
        (function_overall_tension_compression, integration_point_function) = \
            obj_calc.retrieve_functions(selection_tension_compression, integration_point)

        percentile_values = []
        element_count = 0

        for part_id in part_ids:
            index = part_idx[part_id]
            data = part_data[part_id]
            element_ids = np.unique(index[:, 0])
            element_count += len(element_ids)

            for element_id in element_ids:
                row_ids = np.where(index[:, 0] == element_id)[0]

                reduced_integration_point_data = integration_point_function(data[:, row_ids, :], axis=1)
                result_data = function_overall_tension_compression(reduced_integration_point_data, axis=1)
                result_data = result_data[~np.isnan(result_data)]

                # max value for histogram
                max_value_time_index = np.nanargminax(result_data)
                max_value = result_data[max_value_time_index]
                percentile_values.append(max_value)

        # percentile calculation
        percentile_values = sorted(percentile_values)
        if selection_tension_compression == UniversalLimit.TENSION_COMPRESSION_VALUES[2]:
            # compression, flip array to descending order
            percentile_values.reverse()

        value_size = len(percentile_values)
        percentile_value = percentile

        percentile_upper_limit = 1 - (1 / element_count)
        if percentile_value <= percentile_upper_limit:
            percentile_limit = percentile_values[int(np.ceil(percentile_value * value_size))]
        else:
            percentile_limit = percentile_values[int(np.ceil(percentile_upper_limit * value_size))]

        return percentile_limit    

    @staticmethod
    def erf_rib_risk_age(x: pint.Quantity, age: float):
        #x = param_dict['x']
        #age = param_dict['age']
        # parameters from Larsson et al. 2021:
        # https://doi.org/10.3389/fbioe.2021.677768 
        a = 0.3026
        b = -2.9866
        c = -0.013
        return (0.5+0.5*special.erf((np.log((x))-(b+c*age))/(math.sqrt(2)*a)))
        
    @staticmethod
    def extended_ccdf_weibull_ribs_age(x: float, a: float, b: float, c: float, age: float):
        return 1.0 - np.exp(-np.power((x / np.exp(b+c*age)),a))
        
    
    def percentile_part(object_data: dict, selection_tension_compression: str,
                   integration_point: str, percentile: float):
        """only one value possible to return

        Args:
          object_data: param selection_tension_compression:
          integration_point: param percentile:
          object_data: dict: 
          selection_tension_compression: str: 
          integration_point: str: 
          percentile: str: 
          units: 

        Returns:

        """

        part_ids = object_data["part_ids"]
        part_data = object_data["part_data"]
        part_idx = object_data["part_idx"]

        if len(part_ids) == 0:
            return np.NaN

        obj_calc = ObjectCalcUtil()
        (function_overall_tension_compression, integration_point_function) = \
            obj_calc.retrieve_functions(selection_tension_compression, integration_point)


        part_percentile_element = []

        for part_id in part_ids:
            percentile_values = []
            index = part_idx[part_id]
            data = part_data[part_id]
            element_ids = np.unique(index[:, 0])
            element_count = len(element_ids)

            for element_id in element_ids:
                row_ids = np.where(index[:, 0] == element_id)[0]

                reduced_integration_point_data = integration_point_function(data[:, row_ids, :], axis=1)
                result_data = function_overall_tension_compression(reduced_integration_point_data, axis=1)
                result_data = result_data[~np.isnan(result_data)]

                # max value for histogram
                max_value_time_index = np.nanargmax(result_data)
                max_value = result_data[max_value_time_index]
                percentile_values.append(max_value)

            # percentile calculation
            percentile_values = sorted(percentile_values)
            if selection_tension_compression == UniversalLimit.TENSION_COMPRESSION_VALUES[2]:
                # compression, flip array to descending order
                percentile_values.reverse()

            value_size = len(percentile_values)
            percentile_value = percentile

            percentile_upper_limit = 1 - (1 / element_count)
            if percentile_value <= percentile_upper_limit:
                percentile_limit = percentile_values[int(np.ceil(percentile_value * value_size))]
            else:
                percentile_limit = percentile_values[int(np.ceil(percentile_upper_limit * value_size))]
                
            part_percentile_element.append(percentile_limit)

        return part_percentile_element

    @staticmethod
    def OLR(x: pint.Quantity, b0: float, b1: float):
        # ordinal logistic regression
        return 1.0 - (1.0/(1 + np.exp(-(b0 + b1*x))))

    @staticmethod
    def stress_strain_part_min(object_data: dict, selection_tension_compression: str, integration_point: str,
                               nr_smallest_elements: float):
        """

        Args:
          object_data: param selection_tension_compression:
          integration_point: return:
          object_data: dict: 
          selection_tension_compression: str: 
          integration_point: str: 
          nr_smallest_elements: float: 

        Returns:

        """
        part_ids = object_data["part_ids"]
        part_data = object_data["part_data"]
        part_idx = object_data["part_idx"]

        obj_calc = ObjectCalcUtil()
        (function_overall_tension_compression, integration_point_function) = \
            obj_calc.retrieve_functions(selection_tension_compression, integration_point)
        part_min_element = {}

        for part_id in part_ids:
            min_elem = np.empty(shape=(0, 2))

            index = part_idx[part_id]
            data = part_data[part_id]
            if len(index) == 0:
                part_min_element[part_id] = np.empty(int(nr_smallest_elements))
                part_min_element[part_id].fill(np.NaN)
                continue

            element_ids = np.unique(index[:, 0])

            for i, element_id in enumerate(element_ids):
                row_ids = np.where(index[:, 0] == element_id)[0]
                reduced_integration_point_data = integration_point_function(data[:, row_ids, :], axis=1)
                result_data = function_overall_tension_compression(reduced_integration_point_data, axis=1)
                result_data = result_data[~np.isnan(result_data)]

                min_value_time_index = np.nanargmin(result_data)
                min_value = result_data[min_value_time_index]
                min_elem = np.concatenate((min_elem, [[min_value, element_id]]), axis=0)

            ind_of_smallest_elements = min_elem[:, 0].argsort()[-nr_smallest_elements:][::-1]
            min_index = ind_of_smallest_elements[-1]

            part_min_elem_value = min_elem[min_index][0]
            if np.isnan(part_min_elem_value):
                print(LOGConstants.WARNING[0],
                      "Due to deleted elements, the minimum value for " + str(part_id) + " is NaN!")

            part_min_element[part_id] = min_elem[min_index]
            return_array = np.array([part_min_element[i][0] for i in part_min_element])

        return return_array

    @staticmethod
    def HIC_AIS2plus(x: pint.Quantity):
        mean=6.96352
        standard_deviation=0.84664
        return norm.cdf(np.log(x), loc=mean, scale=standard_deviation)

    @staticmethod
    def HIC_AIS3plus(x):
        mean=0
        standard_deviation=1
        return norm.cdf((np.log(x)-7.45231)/0.73998, loc=mean, scale=standard_deviation)

    @staticmethod
    def BrIC_risk(x: pint.Quantity,AIS_factor):
        # based on Takhounts et al. 2013
        # Table III.6 (BrIC based on MPS and formulation 
        # given by equation 4 (average critical angular 
        # velocities from Table 3).
        return 1-np.exp(-np.power((x/AIS_factor),2.84))

    @staticmethod
    def skull_frac_risk(x: pint.Quantity):
        # based on Vander Vorst et al. 2003
        # PMID: 12941236
        e_term = np.exp(5.89*np.log(x)+12.582)
        return e_term/(1+e_term)

    @staticmethod
    def lumbar_vertebra_fracture(force_z: pint.Quantity, moment_y: pint.Quantity, CSA: pint.Quantity, age):
        # based on Östling et al. 2022
        # doi: 10.3389/ffutr.2022.890117
        # expects 
        # 1) a force array (force_z) where compression is positive
        # 2) a moment array (moment_y) where flexion is positive
        # 3) the cross sectional area (CSA) of the main body of the vertebra 
        #    (excluding the posterior region) which is assessed
        # 4) the desired age (age) for which the risk should be calculated
        # if the force/moment arrays do not have the expected polarity, they can be adjusted here
        L = np.zeros(shape=(len(force_z)))
        P = L
        #force_z = force_z * (-1)
        moment_y = moment_y * (-1)
        for i in range(len(L)):
            L = (1 - 0.11)*force_z[i].to(ureg.parse_expression('newton')).magnitude/(CSA.to(ureg.parse_expression('mm**2')).magnitude) \
                + 0.11*moment_y[i].to(ureg.parse_expression('newton*millimeter')).magnitude/np.power((CSA.to(ureg.parse_expression('mm**2')).magnitude),3/2)
            if math.isnan(1-np.exp(-np.power(L/(np.exp(1.896-0.00874*age)),(1/0.201)))):
                P[i]=0
            else:
                P[i]=1-np.exp(-np.power(L/(np.exp(1.896-0.00874*age)),(1/0.201)))
        return P

    @staticmethod
    def defl_based_NFR3plus_risk(r_max: pint.Quantity, age):
        # based on Craig et al. 2020
        # function below expects maximum deflection to be a positive number
        # hence, the absolute value of the negative deflection is used - adjust if necessary.
        r_max = abs(r_max).to(ureg.parse_expression('mm')).magnitude
        p=1-np.exp(-((r_max/np.exp(4.7276-0.0166*age))**2.977))
        return p

    @staticmethod
    def defl_based_NFR7plus_risk(r_max: pint.Quantity, age):
        # based on Craig et al. 2020
        # function below expects maximum deflection to be a positive number
        # hence, the absolute value of the negative deflection is used - adjust if necessary.
        r_max = abs(r_max).to(ureg.parse_expression('mm')).magnitude
        p=1-np.exp(-((r_max/np.exp(4.4272-0.0077*age))**3.9245))
        return p

# use to define own functions class UserFunction: @staticmethod def time_at_Xmove(param_dict, units): move = (param_dict['movement']) time = (param_dict['time']) value = (param_dict['value']) offset = (param_dict['offset']) timestep_of_maximum = 0 for t in range(len(move)): y = move[t] if y > value: timestep_of_maximum = t break time_of_maximum = time[timestep_of_maximum] - offset return (np.max(time_of_maximum)) @staticmethod def vel_at_timeX(param_dict, units): move = (param_dict['movement']) vel = (param_dict['velocity']) value = (param_dict['value']) for t in range(len(move)): y = move[t] if y > value: timestep_of_maximum = t break velocity_at_timeX = vel[timestep_of_maximum] return (np.max(velocity_at_timeX)) @staticmethod def FC_t(param_dict, units): t = param_dict['time'].flatten() # / units.second() # t in [s] fr = param_dict['Fres'].flatten() fbmax = max(fr) f_rise = fbmax * 0.01 trise = -1 for i in range(len(t) - 1): if f_rise < fr[i]: trise = t[i] break # inclusive 20ms start_time! -- zeitpunkt in der fbelt um 20% gefallen ist return trise @staticmethod def H350_FEMUR_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 # FEMUR scalar value - injuryvalue risk = 1 / (1 + math.exp(offset - scaling * (np.abs(injuryvalue)))) return (risk) @staticmethod def H350_NECKCOMP_AIS3(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * (np.abs(injuryvalue[0])))) return (risk) @staticmethod def H350_NECKTEN_AIS3(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * (np.abs(injuryvalue[0])))) return (risk) @staticmethod def H350_CHESTDEF_AIS3(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * ((np.abs(injuryvalue[0])) ** power))) return (risk) @staticmethod def H350_NIJ_AIS3(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * injuryvalue[0])) return (risk) @staticmethod def THOR50_NIJ_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * injuryvalue[0])) return (risk) @staticmethod def THOR50_NIJ_AIS3(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * injuryvalue[0])) # "name" : "THOR50_NIJ_AIS3+", # "ccdf" : "1 / (1 + math.exp(4.9372 - 4.5294*value))" return (risk) @staticmethod def THOR50_CHESTDEF_AIS3(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 - math.exp(-((injuryvalue / divide) ** power)) # "name" : "THOR50_CHESTDEF_AIS3+", # "ccdf" : "1 - math.exp(-((value/59.865)**2.7187))" return (risk) @staticmethod def THOR50_ABDOMENDEF_AIS3(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * injuryvalue)) # "name" : "THOR50_ABDOMENDEF_AIS3+", # "ccdf" : "1 / (1 + math.exp(7.849 - 0.0886 * value))" return (risk) @staticmethod def THOR50_FEMUR_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * (abs(injuryvalue)))) # "name" : "THOR50_FEMUR_AIS2+", # "ccdf" : "1 / (1 + math.exp(5.7949 - 0.6748*(abs(value))))" return (risk) @staticmethod def THOR50_TIBIA_FZ_UPPER_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * injuryvalue)) # "name" : "THOR50_TIBIA_FZ_UPPER_AIS2+", # "ccdf" : "1 / (1 + math.exp(5.6654 - 0.8189*value))" return (risk) @staticmethod def THOR50_TIBIA_FZ_LOWER_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 / (1 + math.exp(offset - scaling * injuryvalue)) # "name" : "THOR50_TIBIA_FZ_LOWER_AIS2+", # "ccdf" : "1 / (1 + math.exp(3.9121 - 0.48*value))" return (risk) @staticmethod def THOR50_TIBIA_Mres_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = 1 - math.exp((math.log(injuryvalue) - offset) / divide) # "name" : "THOR50_TIBIA_Mres_AIS2+", # "ccdf" : "1 - math.exp( - math.exp((math.log(value)-5.8492)/0.2965))" return (risk) @staticmethod def THOR50_HIP_Fracture_AIS(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] risk = 0 risk = (1.0 + math.erf(((math.log(scaling * injuryvalue[0]) - offset) / divide) / math.sqrt(2.0))) / 2.0 # "name" : "THOR50_HIP_Fracture_AIS+", # "ccdf" : "(1.0 + math.erf(((math.log(1.429*value) - 1.6058) / 0.2339) / math.sqrt(2.0))) / 2.0" return (risk) @staticmethod def WS50_SHOULDERFORCE_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] age = param_dict["age"] risk = 0 risk = 1 - math.exp(-(np.abs(injuryvalue) / (offset - scaling * age)) ** power) # "name" : "WS50_SHOULDERFORCE_AIS2+", # "ccdf" : "1 - (math.exp(-(value / (8.144 - 0.006 * 50)) ** 7.41))" return (risk) @staticmethod def WS50_ABDOMENDEF_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] age = param_dict["age"] risk = 0 risk = 1 - math.exp(-((np.abs(injuryvalue[0]) / (offset - scaling * age))) ** power) # "name" : "WS50_ABDOMENDEF_AIS2+" Soft Tissue Abdominal Injury, # "ccdf" : "1 - math.exp( - (injuryvalue[0] / (5.368 - 0.021 * 50)) ** 8.61)" return (risk) @staticmethod def WS50_DEF_SKEL_AIS3(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] age = param_dict["age"] risk = 0 risk = 1 / (1 + math.exp( -(math.log(abs(injuryvalue[0])) - (offset - scaling * age)) / divide)) return (risk) @staticmethod def WS50_PUBICFORCE_AIS2(param_dict, units): injuryvalue = param_dict["injuryvalue"] offset = param_dict["offset"] divide = param_dict["divide"] scaling = param_dict["scaling"] power = param_dict["power"] age = param_dict["age"] risk = 0 risk = 1 - math.exp( - ( ( injuryvalue[0] / (offset - scaling * age) ) ** power) ) # "name" : "WS50_PUBICFORCE_AIS2+", # "ccdf" : "1 - math.exp(-(value / (8.775 - 0.014 * 50)) ** 4.60)" return (risk) @staticmethod def THOR_TCI_Rmax(Rmax_Points_UL: pint.Quantity, Rmax_Points_LL: pint.Quantity, Rmax_Points_UR: pint.Quantity, Rmax_Points_LR: pint.Quantity): Rmax = 0 * Rmax_Points_UL.units Rmax = max([np.max(Rmax_Points_UL[:, 0]), np.max(Rmax_Points_LL[:, 0]), np.max(Rmax_Points_UR[:, 0]), np.max(Rmax_Points_LR[:, 0])]) return (Rmax) @staticmethod def THOR_Abdomen_Deflection(Abdmax_LEFT: pint.Quantity, Abdmax_RIGHT: pint.Quantity): Abdomenmax = 0 * Abdmax_LEFT.units Abdomenmax = max([np.max(Abdmax_RIGHT[:, 0]), np.max(Abdmax_LEFT[:, 0])]) return (Abdomenmax) @staticmethod def AIRBAG_CONTACT(param_dict, units): data_vector = param_dict['data_vector'] time = param_dict['time'] move = param_dict['move'] tbase = param_dict['tbase'] t_cont = 0 for i in range(len(time) - 1): if move < np.abs(data_vector[i]): t_cont = time[i] - tbase break return np.max(t_cont) @staticmethod def VEL_NODE_MAX(param_dict, units): t = param_dict['time'].flatten() y = param_dict['y'].flatten() vel = [] velmax = 0 vel.append(velmax) for i in range(len(t) - 1): if i > 0: vel.append((y[i] - y[i - 1]) / (t[i] - t[i - 1])) if np.abs(vel[i]) > np.abs(velmax): velmax = vel[i] return (velmax) @staticmethod def VEL_NODE(param_dict, units): t = param_dict['time'].flatten() y = param_dict['y'].flatten() vel = [] velmax = 0 vel.append(velmax) for i in range(len(t) - 1): if i > 0: vel.append((y[i] - y[i - 1]) / (t[i] - t[i - 1])) return np.array(vel).reshape(-1, 1) @staticmethod def ext_bend_mom(param_dict, units): force_x = param_dict['force_x'] / (units.meter() / (units.second() ** 2)) moment_y = param_dict['moment_y'] / ((units.meter() ** 2) / (units.second() ** 2)) moc_d = param_dict['distance_occipital_condyle'] / ((units.meter() ** 2) / (units.second() ** 2)) mtot = moment_y - moc_d * force_x return (mtot) @staticmethod def double_abs_sub_res(arr1: pint.Quantity, arr2: pint.Quantity, arr3: pint.Quantity, arr4: float): arr1m = arr1.magnitude arr2m = arr2.magnitude arr3m = arr3.magnitude arr4m = arr4.magnitude length = arr3.units rel_dist_max_1 = [] rel_dist_max_2 = [] rel_dist = [] rel_dist_1 = np.subtract(arr1m, arr2m) rel_dist_2 = np.subtract(arr1m, arr3m) for i in range(len(rel_dist_1)): rel_dist_max_1.append( (rel_dist_1[i, 0] ** 2 + rel_dist_1[i, 1] ** 2 + rel_dist_1[i, 2] ** 2) ** (1 / 2) - arr4m) for i in range(len(rel_dist_2)): rel_dist_max_2.append( (rel_dist_2[i, 0] ** 2 + rel_dist_2[i, 1] ** 2 + rel_dist_2[i, 2] ** 2) ** (1 / 2) - arr4m) if np.min(rel_dist_max_1) < 0 or np.min(rel_dist_max_2) < 0: for i in range(len(rel_dist_1)): rel_dist.append(0) else: if np.max(np.abs(rel_dist_max_1)) < np.max(np.abs(rel_dist_max_2)): rel_dist = rel_dist_max_1 else: rel_dist = rel_dist_max_2 return np.abs(rel_dist) * length @staticmethod def MOCy_max(param_dict, units): force_x = param_dict['force_x'] / (units.meter() / (units.second() ** 2)) moment_y = param_dict['moment_y'] / ((units.meter() ** 2) / (units.second() ** 2)) moc_d = param_dict['distance_occipital_condyle'] MOCy_max = np.max(moment_y - moc_d * force_x) return (MOCy_max) @staticmethod def MOCy_min(param_dict, units): force_x = param_dict['force_x'] / (units.meter() / (units.second() ** 2)) moment_y = param_dict['moment_y'] / ((units.meter() ** 2) / (units.second() ** 2)) moc_d = param_dict['distance_occipital_condyle'] MOCy_min = np.min(moment_y - moc_d * force_x) return (MOCy_min) @staticmethod def minimum_of_four(data_vector_1: pint.Quantity, data_vector_2: pint.Quantity, data_vector_3: pint.Quantity, data_vector_4: pint.Quantity): min_1 = (np.min(data_vector_1)).flatten() min_2 = (np.min(data_vector_2)).flatten() min_3 = (np.min(data_vector_3)).flatten() min_4 = (np.min(data_vector_4)).flatten() temp = [min_1, min_2, min_3, min_4] min = np.min(temp) # index_of_minimum = np.argmin(param_dict['minimum']) # time_of_minimum = param_dict['time'][index_of_minimum] return (min) * data_vector_1.units @staticmethod def maximum_of_four(data_vector_1: pint.Quantity, data_vector_2: pint.Quantity, data_vector_3: pint.Quantity, data_vector_4: pint.Quantity): max_1 = (np.max(data_vector_1)).flatten() max_2 = (np.max(data_vector_2)).flatten() max_3 = (np.max(data_vector_3)).flatten() max_4 = (np.max(data_vector_4)).flatten() temp = [max_1, max_2, max_3, max_4] max = np.max(temp) # index_of_minimum = np.argmin(param_dict['minimum']) # time_of_minimum = param_dict['time'][index_of_minimum] return (max) * data_vector_1.units @staticmethod def ABS_MAX_VALUE(param_dict, units): data_vec = param_dict['data_vector'] val = 0 if np.abs(np.min(data_vec)) < np.abs(np.max(data_vec)): val = np.max(data_vec) else: val = np.min(data_vec) return val @staticmethod def RES_ACC_from_XYZ_Velocities(a_x: pint.Quantity, a_y: pint.Quantity, a_z: pint.Quantity): a_x = param_dict['a_x'].flatten() a_y = param_dict['a_y'].flatten() a_z = param_dict['a_z'].flatten() a_res = (a_x ** 2 + a_y ** 2 + a_z ** 2) ** (1 / 2) return np.array(a_res).reshape(-1, 1) @staticmethod def ILIAC_Drop(Fz: pint.Quantity, time: pint.Quantity): t = time.to_base_units().magnitude fi = Fz.magnitude iliac_drop = [] delta_t = (t[1] - t[0])[0] tresample = int(round(0.001 / delta_t)) for i in range(len(t) - (1 + tresample)): value = fi[i + tresample] - fi[i] if delta_t < 0.001 and value < 0: iliac_drop.append(value[0]) else: iliac_drop.append(0) return np.array(iliac_drop).reshape(-1, 1) * time.units @staticmethod def FB_drop_t(time: pint.Quantity, Fb: pint.Quantity, drop_value_percent: pint.Quantity): t = time fb = Fb dvp = drop_value_percent fbmax = max(fb) fbmax_drop = fbmax * (1 - dvp / 100) tswitch = 0 fb_avg = 0 for i in range(len(t) - 1): fb_avg = fb_avg + fb[i] if fbmax == fb[i]: tmax = t[i] tswitch = 1 if tswitch == 1: if fbmax_drop >= fb[i]: tdrop = t[i] average_i = i tswitch = 0 break if tswitch == 0: tdrop = t[i] average_i = i # durchschnittliche Kraft bis fbmax_drop fb_avg_drop = fb_avg / average_i # inclusive 20ms start_time! -- zeitpunkt in der fbelt um 20% gefallen ist return tdrop @staticmethod def FB_drop_avg(Fb: pint.Quantity, time: pint.Quantity, drop_value_percent: pint.Quantity): t = time fb = Fb dvp = drop_value_percent.magnitude fbmax = max(fb) fbmax_drop = fbmax * (1 - dvp / 100) tswitch = 0 fb_avg = 0 for i in range(len(t) - 1): fb_avg = fb_avg + fb[i] if fbmax == fb[i]: tmax = t[i] tswitch = 1 if tswitch == 1: if fbmax_drop >= fb[i]: tdrop = t[i] average_i = i tswitch = 0 break if tswitch == 0: tdrop = t[i] average_i = i # durchschnittliche Kraft bis fbmax_drop fb_avg_drop = fb_avg / average_i # nur ein return return fb_avg_drop # END Peter L. - BeltLoad80%########### @staticmethod def HIC_AIS2(HIC: pint.Quantity): HICAIS2 = norm.cdf((math.log(HIC) - 6.96352) / 0.84664) return HICAIS2 @staticmethod def HIC_AIS3(HIC: pint.Quantity): HICAIS3 = norm.cdf((math.log(HIC) - 7.45231) / 0.73998) return HICAIS3 @staticmethod def HIC_AIS4(HIC: pint.Quantity): HICAIS4 = norm.cdf((math.log(HIC) - 7.65605) / 0.60580) return HICAIS4 @staticmethod def BRIC_AIS4(BRIC: pint.Quantity): BRICAIS4 = 1 - np.exp((BRIC-0.523)/0.647) ** 1.8 return BRICAIS4 @staticmethod def NIJ_AIS2_THOR(NIJ: pint.Quantity): NIJAIS2 = 1 / (1 + np.exp(4.3085 - 5.4079 * NIJ)) return NIJAIS2 @staticmethod def NIJ_AIS3_THOR(NIJ: pint.Quantity): NIJAIS3 = 1 / (1 + np.exp(4.9372 - 4.5294 * NIJ)) return NIJAIS3 @staticmethod def NIJ_AIS3_H305(NIJ: pint.Quantity): NIJAIS3 = 1 - np.exp(-(NIJ/1.3933) ** 2.8816) return NIJAIS3 @staticmethod def THOR_TCI_Rmax_AIS3(TCI_Rmax: pint.Quantity): TCI_RmaxAIS3 = 1 - np.exp(-(THOR_TCI_Rmax_AIS3/59.865) ** 2.7187) return TCI_RmaxAIS3 @staticmethod def Abdomen_Deflection_AIS3_THOR(Abdomen_Deflection: pint.Quantity): Abdomen_DeflectionAIS3 = 1 / (1 + np.exp(7.849 - 0.0886 * Abdomen_Deflection)) return Abdomen_DeflectionAIS3 @staticmethod def Femur_axial_force_AIS2_THOR(Femur_Comp: pint.Quantity): Femur_Comp_mag = Femur_Comp.magnitude Femur_axialAIS2 = 1 / (1 + np.exp(5.7949 - 0.6748 * Femur_Comp_mag)) return Femur_axialAIS2 @staticmethod def Femur_axial_force_AIS2_H305(Femur_Comp: pint.Quantity): Femur_Comp_mag = Femur_Comp.magnitude Femur_axialAIS2 = 1 / (1 + np.exp(5.7949 - 0.7619 * Femur_Comp_mag)) return Femur_axialAIS2 @staticmethod def Lower_Tibia_axial_force_AIS2_THOR(Tibia_Fz: pint.Quantity): Tibia_Fz_mag = Tibia_Fz.magnitude Tibia_FzAIS2 = 1 / (1 + np.exp(3.9121 - 0.48 * Tibia_Fz_mag)) return Tibia_FzAIS2 @staticmethod def Upper_Tibia_axial_force_AIS2_THOR(Tibia_Fz: pint.Quantity): Tibia_Fz_mag = Tibia_Fz.magnitude Tibia_FzAIS2 = 1 / (1 + np.exp(5.6654 - 0.8189 * Tibia_Fz_mag)) return Tibia_FzAIS2 @staticmethod def Tibia_moment_AIS2_THOR(Tibia_Moment: pint.Quantity): Tibia_Momentm =Tibia_Moment.magnitude Tibia_MomentAIS2 = 1 - np.exp(-np.exp((math.log(Tibia_Momentm) - 5.8492)/0.2965)) return Tibia_MomentAIS2 @staticmethod def initValue (arr1: pint.Quantity): initialValue = arr1[0] return initialValue @staticmethod def res_wo_unit(arr1: pint.Quantity, arr2: pint.Quantity, arr3: pint.Quantity): arr1m = arr1.magnitude arr2m = arr2.magnitude arr3m = arr3.magnitude unit = arr1.units return np.sqrt(arr1m**2 + arr2m**2 + arr3m**2) * unit @staticmethod def percentilemin(object_data: dict, selection_tension_compression: str, integration_point: str, percentile: float): """only one value possible to return Args: object_data: param selection_tension_compression: integration_point: param percentile: object_data: dict: selection_tension_compression: str: integration_point: str: percentile: str: units: Returns: """ part_ids = object_data["part_ids"] part_data = object_data["part_data"] part_idx = object_data["part_idx"] if len(part_ids) == 0: return np.NaN obj_calc = ObjectCalcUtil() (function_overall_tension_compression, integration_point_function) = \ obj_calc.retrieve_functions(selection_tension_compression, integration_point) percentile_values = [] element_count = 0 for part_id in part_ids: index = part_idx[part_id] data = part_data[part_id] element_ids = np.unique(index[:, 0]) element_count += len(element_ids) for element_id in element_ids: row_ids = np.where(index[:, 0] == element_id)[0] reduced_integration_point_data = integration_point_function(data[:, row_ids, :], axis=1) result_data = function_overall_tension_compression(reduced_integration_point_data, axis=1) result_data = result_data[~np.isnan(result_data)] # max value for histogram max_value_time_index = np.nanargminax(result_data) max_value = result_data[max_value_time_index] percentile_values.append(max_value) # percentile calculation percentile_values = sorted(percentile_values) if selection_tension_compression == UniversalLimit.TENSION_COMPRESSION_VALUES[2]: # compression, flip array to descending order percentile_values.reverse() value_size = len(percentile_values) percentile_value = percentile percentile_upper_limit = 1 - (1 / element_count) if percentile_value <= percentile_upper_limit: percentile_limit = percentile_values[int(np.ceil(percentile_value * value_size))] else: percentile_limit = percentile_values[int(np.ceil(percentile_upper_limit * value_size))] return percentile_limit @staticmethod def erf_rib_risk_age(x: pint.Quantity, age: float): #x = param_dict['x'] #age = param_dict['age'] # parameters from Larsson et al. 2021: # https://doi.org/10.3389/fbioe.2021.677768 a = 0.3026 b = -2.9866 c = -0.013 return (0.5+0.5*special.erf((np.log((x))-(b+c*age))/(math.sqrt(2)*a))) @staticmethod def extended_ccdf_weibull_ribs_age(x: float, a: float, b: float, c: float, age: float): return 1.0 - np.exp(-np.power((x / np.exp(b+c*age)),a)) def percentile_part(object_data: dict, selection_tension_compression: str, integration_point: str, percentile: float): """only one value possible to return Args: object_data: param selection_tension_compression: integration_point: param percentile: object_data: dict: selection_tension_compression: str: integration_point: str: percentile: str: units: Returns: """ part_ids = object_data["part_ids"] part_data = object_data["part_data"] part_idx = object_data["part_idx"] if len(part_ids) == 0: return np.NaN obj_calc = ObjectCalcUtil() (function_overall_tension_compression, integration_point_function) = \ obj_calc.retrieve_functions(selection_tension_compression, integration_point) part_percentile_element = [] for part_id in part_ids: percentile_values = [] index = part_idx[part_id] data = part_data[part_id] element_ids = np.unique(index[:, 0]) element_count = len(element_ids) for element_id in element_ids: row_ids = np.where(index[:, 0] == element_id)[0] reduced_integration_point_data = integration_point_function(data[:, row_ids, :], axis=1) result_data = function_overall_tension_compression(reduced_integration_point_data, axis=1) result_data = result_data[~np.isnan(result_data)] # max value for histogram max_value_time_index = np.nanargmax(result_data) max_value = result_data[max_value_time_index] percentile_values.append(max_value) # percentile calculation percentile_values = sorted(percentile_values) if selection_tension_compression == UniversalLimit.TENSION_COMPRESSION_VALUES[2]: # compression, flip array to descending order percentile_values.reverse() value_size = len(percentile_values) percentile_value = percentile percentile_upper_limit = 1 - (1 / element_count) if percentile_value <= percentile_upper_limit: percentile_limit = percentile_values[int(np.ceil(percentile_value * value_size))] else: percentile_limit = percentile_values[int(np.ceil(percentile_upper_limit * value_size))] part_percentile_element.append(percentile_limit) return part_percentile_element @staticmethod def OLR(x: pint.Quantity, b0: float, b1: float): # ordinal logistic regression return 1.0 - (1.0/(1 + np.exp(-(b0 + b1*x)))) @staticmethod def stress_strain_part_min(object_data: dict, selection_tension_compression: str, integration_point: str, nr_smallest_elements: float): """ Args: object_data: param selection_tension_compression: integration_point: return: object_data: dict: selection_tension_compression: str: integration_point: str: nr_smallest_elements: float: Returns: """ part_ids = object_data["part_ids"] part_data = object_data["part_data"] part_idx = object_data["part_idx"] obj_calc = ObjectCalcUtil() (function_overall_tension_compression, integration_point_function) = \ obj_calc.retrieve_functions(selection_tension_compression, integration_point) part_min_element = {} for part_id in part_ids: min_elem = np.empty(shape=(0, 2)) index = part_idx[part_id] data = part_data[part_id] if len(index) == 0: part_min_element[part_id] = np.empty(int(nr_smallest_elements)) part_min_element[part_id].fill(np.NaN) continue element_ids = np.unique(index[:, 0]) for i, element_id in enumerate(element_ids): row_ids = np.where(index[:, 0] == element_id)[0] reduced_integration_point_data = integration_point_function(data[:, row_ids, :], axis=1) result_data = function_overall_tension_compression(reduced_integration_point_data, axis=1) result_data = result_data[~np.isnan(result_data)] min_value_time_index = np.nanargmin(result_data) min_value = result_data[min_value_time_index] min_elem = np.concatenate((min_elem, [[min_value, element_id]]), axis=0) ind_of_smallest_elements = min_elem[:, 0].argsort()[-nr_smallest_elements:][::-1] min_index = ind_of_smallest_elements[-1] part_min_elem_value = min_elem[min_index][0] if np.isnan(part_min_elem_value): print(LOGConstants.WARNING[0], "Due to deleted elements, the minimum value for " + str(part_id) + " is NaN!") part_min_element[part_id] = min_elem[min_index] return_array = np.array([part_min_element[i][0] for i in part_min_element]) return return_array @staticmethod def HIC_AIS2plus(x: pint.Quantity): mean=6.96352 standard_deviation=0.84664 return norm.cdf(np.log(x), loc=mean, scale=standard_deviation) @staticmethod def HIC_AIS3plus(x): mean=0 standard_deviation=1 return norm.cdf((np.log(x)-7.45231)/0.73998, loc=mean, scale=standard_deviation) @staticmethod def BrIC_risk(x: pint.Quantity,AIS_factor): # based on Takhounts et al. 2013 # Table III.6 (BrIC based on MPS and formulation # given by equation 4 (average critical angular # velocities from Table 3). return 1-np.exp(-np.power((x/AIS_factor),2.84)) @staticmethod def skull_frac_risk(x: pint.Quantity): # based on Vander Vorst et al. 2003 # PMID: 12941236 e_term = np.exp(5.89*np.log(x)+12.582) return e_term/(1+e_term) @staticmethod def lumbar_vertebra_fracture(force_z: pint.Quantity, moment_y: pint.Quantity, CSA: pint.Quantity, age): # based on Östling et al. 2022 # doi: 10.3389/ffutr.2022.890117 # expects # 1) a force array (force_z) where compression is positive # 2) a moment array (moment_y) where flexion is positive # 3) the cross sectional area (CSA) of the main body of the vertebra # (excluding the posterior region) which is assessed # 4) the desired age (age) for which the risk should be calculated # if the force/moment arrays do not have the expected polarity, they can be adjusted here L = np.zeros(shape=(len(force_z))) P = L #force_z = force_z * (-1) moment_y = moment_y * (-1) for i in range(len(L)): L = (1 - 0.11)*force_z[i].to(ureg.parse_expression('newton')).magnitude/(CSA.to(ureg.parse_expression('mm**2')).magnitude) \ + 0.11*moment_y[i].to(ureg.parse_expression('newton*millimeter')).magnitude/np.power((CSA.to(ureg.parse_expression('mm**2')).magnitude),3/2) if math.isnan(1-np.exp(-np.power(L/(np.exp(1.896-0.00874*age)),(1/0.201)))): P[i]=0 else: P[i]=1-np.exp(-np.power(L/(np.exp(1.896-0.00874*age)),(1/0.201))) return P @staticmethod def defl_based_NFR3plus_risk(r_max: pint.Quantity, age): # based on Craig et al. 2020 # function below expects maximum deflection to be a positive number # hence, the absolute value of the negative deflection is used - adjust if necessary. r_max = abs(r_max).to(ureg.parse_expression('mm')).magnitude p=1-np.exp(-((r_max/np.exp(4.7276-0.0166*age))**2.977)) return p @staticmethod def defl_based_NFR7plus_risk(r_max: pint.Quantity, age): # based on Craig et al. 2020 # function below expects maximum deflection to be a positive number # hence, the absolute value of the negative deflection is used - adjust if necessary. r_max = abs(r_max).to(ureg.parse_expression('mm')).magnitude p=1-np.exp(-((r_max/np.exp(4.4272-0.0077*age))**3.9245)) return p