text:
#- Other VISA Definitions --------------------------------------------------*/
const WORD_SIZE = 64
const VI_NULL = 0

const VI_TRUE = 1
const VI_FALSE = 0

#- Attributes (platform independent size) ----------------------------------*/

const VI_ATTR_RSRC_CLASS                  =  0xBFFF0001
const VI_ATTR_RSRC_NAME                   =  0xBFFF0002
const VI_ATTR_RSRC_IMPL_VERSION           =  0x3FFF0003
const VI_ATTR_RSRC_LOCK_STATE             =  0x3FFF0004
const VI_ATTR_MAX_QUEUE_LENGTH            =  0x3FFF0005
const VI_ATTR_USER_DATA_32                =  0x3FFF0007
const VI_ATTR_FDC_CHNL                    =  0x3FFF000D
const VI_ATTR_FDC_MODE                    =  0x3FFF000F
const VI_ATTR_FDC_GEN_SIGNAL_EN           =  0x3FFF0011
const VI_ATTR_FDC_USE_PAIR                =  0x3FFF0013
const VI_ATTR_SEND_END_EN                 =  0x3FFF0016
const VI_ATTR_TERMCHAR                    =  0x3FFF0018
const VI_ATTR_TMO_VALUE                   =  0x3FFF001A
const VI_ATTR_GPIB_READDR_EN              =  0x3FFF001B
const VI_ATTR_IO_PROT                     =  0x3FFF001C
const VI_ATTR_DMA_ALLOW_EN                =  0x3FFF001E
const VI_ATTR_ASRL_BAUD                   =  0x3FFF0021
const VI_ATTR_ASRL_DATA_BITS              =  0x3FFF0022
const VI_ATTR_ASRL_PARITY                 =  0x3FFF0023
const VI_ATTR_ASRL_STOP_BITS              =  0x3FFF0024
const VI_ATTR_ASRL_FLOW_CNTRL             =  0x3FFF0025
const VI_ATTR_RD_BUF_OPER_MODE            =  0x3FFF002A
const VI_ATTR_RD_BUF_SIZE                 =  0x3FFF002B
const VI_ATTR_WR_BUF_OPER_MODE            =  0x3FFF002D
const VI_ATTR_WR_BUF_SIZE                 =  0x3FFF002E
const VI_ATTR_SUPPRESS_END_EN             =  0x3FFF0036
const VI_ATTR_TERMCHAR_EN                 =  0x3FFF0038
const VI_ATTR_DEST_ACCESS_PRIV            =  0x3FFF0039
const VI_ATTR_DEST_BYTE_ORDER             =  0x3FFF003A
const VI_ATTR_SRC_ACCESS_PRIV             =  0x3FFF003C
const VI_ATTR_SRC_BYTE_ORDER              =  0x3FFF003D
const VI_ATTR_SRC_INCREMENT               =  0x3FFF0040
const VI_ATTR_DEST_INCREMENT              =  0x3FFF0041
const VI_ATTR_WIN_ACCESS_PRIV             =  0x3FFF0045
const VI_ATTR_WIN_BYTE_ORDER              =  0x3FFF0047
const VI_ATTR_GPIB_ATN_STATE              =  0x3FFF0057
const VI_ATTR_GPIB_ADDR_STATE             =  0x3FFF005C
const VI_ATTR_GPIB_CIC_STATE              =  0x3FFF005E
const VI_ATTR_GPIB_NDAC_STATE             =  0x3FFF0062
const VI_ATTR_GPIB_SRQ_STATE              =  0x3FFF0067
const VI_ATTR_GPIB_SYS_CNTRL_STATE        =  0x3FFF0068
const VI_ATTR_GPIB_HS488_CBL_LEN          =  0x3FFF0069
const VI_ATTR_CMDR_LA                     =  0x3FFF006B
const VI_ATTR_VXI_DEV_CLASS               =  0x3FFF006C
const VI_ATTR_MAINFRAME_LA                =  0x3FFF0070
const VI_ATTR_MANF_NAME                   =  0xBFFF0072
const VI_ATTR_MODEL_NAME                  =  0xBFFF0077
const VI_ATTR_VXI_VME_INTR_STATUS         =  0x3FFF008B
const VI_ATTR_VXI_TRIG_STATUS             =  0x3FFF008D
const VI_ATTR_VXI_VME_SYSFAIL_STATE       =  0x3FFF0094
const VI_ATTR_WIN_BASE_ADDR_32            =  0x3FFF0098
const VI_ATTR_WIN_SIZE_32                 =  0x3FFF009A
const VI_ATTR_ASRL_AVAIL_NUM              =  0x3FFF00AC
const VI_ATTR_MEM_BASE_32                 =  0x3FFF00AD
const VI_ATTR_ASRL_CTS_STATE              =  0x3FFF00AE
const VI_ATTR_ASRL_DCD_STATE              =  0x3FFF00AF
const VI_ATTR_ASRL_DSR_STATE              =  0x3FFF00B1
const VI_ATTR_ASRL_DTR_STATE              =  0x3FFF00B2
const VI_ATTR_ASRL_END_IN                 =  0x3FFF00B3
const VI_ATTR_ASRL_END_OUT                =  0x3FFF00B4
const VI_ATTR_ASRL_REPLACE_CHAR           =  0x3FFF00BE
const VI_ATTR_ASRL_RI_STATE               =  0x3FFF00BF
const VI_ATTR_ASRL_RTS_STATE              =  0x3FFF00C0
const VI_ATTR_ASRL_XON_CHAR               =  0x3FFF00C1
const VI_ATTR_ASRL_XOFF_CHAR              =  0x3FFF00C2
const VI_ATTR_WIN_ACCESS                  =  0x3FFF00C3
const VI_ATTR_RM_SESSION                  =  0x3FFF00C4
const VI_ATTR_VXI_LA                      =  0x3FFF00D5
const VI_ATTR_MANF_ID                     =  0x3FFF00D9
const VI_ATTR_MEM_SIZE_32                 =  0x3FFF00DD
const VI_ATTR_MEM_SPACE                   =  0x3FFF00DE
const VI_ATTR_MODEL_CODE                  =  0x3FFF00DF
const VI_ATTR_SLOT                        =  0x3FFF00E8
const VI_ATTR_INTF_INST_NAME              =  0xBFFF00E9
const VI_ATTR_IMMEDIATE_SERV              =  0x3FFF0100
const VI_ATTR_INTF_PARENT_NUM             =  0x3FFF0101
const VI_ATTR_RSRC_SPEC_VERSION           =  0x3FFF0170
const VI_ATTR_INTF_TYPE                   =  0x3FFF0171
const VI_ATTR_GPIB_PRIMARY_ADDR           =  0x3FFF0172
const VI_ATTR_GPIB_SECONDARY_ADDR         =  0x3FFF0173
const VI_ATTR_RSRC_MANF_NAME              =  0xBFFF0174
const VI_ATTR_RSRC_MANF_ID                =  0x3FFF0175
const VI_ATTR_INTF_NUM                    =  0x3FFF0176
const VI_ATTR_TRIG_ID                     =  0x3FFF0177
const VI_ATTR_GPIB_REN_STATE              =  0x3FFF0181
const VI_ATTR_GPIB_UNADDR_EN              =  0x3FFF0184
const VI_ATTR_DEV_STATUS_BYTE             =  0x3FFF0189
const VI_ATTR_FILE_APPEND_EN              =  0x3FFF0192
const VI_ATTR_VXI_TRIG_SUPPORT            =  0x3FFF0194
const VI_ATTR_TCPIP_ADDR                  =  0xBFFF0195
const VI_ATTR_TCPIP_HOSTNAME              =  0xBFFF0196
const VI_ATTR_TCPIP_PORT                  =  0x3FFF0197
const VI_ATTR_TCPIP_DEVICE_NAME           =  0xBFFF0199
const VI_ATTR_TCPIP_NODELAY               =  0x3FFF019A
const VI_ATTR_TCPIP_KEEPALIVE             =  0x3FFF019B
const VI_ATTR_4882_COMPLIANT              =  0x3FFF019F
const VI_ATTR_USB_SERIAL_NUM              =  0xBFFF01A0
const VI_ATTR_USB_INTFC_NUM               =  0x3FFF01A1
const VI_ATTR_USB_PROTOCOL                =  0x3FFF01A7
const VI_ATTR_USB_MAX_INTR_SIZE           =  0x3FFF01AF
const VI_ATTR_PXI_DEV_NUM                 =  0x3FFF0201
const VI_ATTR_PXI_FUNC_NUM                =  0x3FFF0202
const VI_ATTR_PXI_BUS_NUM                 =  0x3FFF0205
const VI_ATTR_PXI_CHASSIS                 =  0x3FFF0206
const VI_ATTR_PXI_SLOTPATH                =  0xBFFF0207
const VI_ATTR_PXI_SLOT_LBUS_LEFT          =  0x3FFF0208
const VI_ATTR_PXI_SLOT_LBUS_RIGHT         =  0x3FFF0209
const VI_ATTR_PXI_TRIG_BUS                =  0x3FFF020A
const VI_ATTR_PXI_STAR_TRIG_BUS           =  0x3FFF020B
const VI_ATTR_PXI_STAR_TRIG_LINE          =  0x3FFF020C
const VI_ATTR_PXI_SRC_TRIG_BUS            =  0x3FFF020D
const VI_ATTR_PXI_DEST_TRIG_BUS           =  0x3FFF020E
const VI_ATTR_PXI_MEM_TYPE_BAR0           =  0x3FFF0211
const VI_ATTR_PXI_MEM_TYPE_BAR1           =  0x3FFF0212
const VI_ATTR_PXI_MEM_TYPE_BAR2           =  0x3FFF0213
const VI_ATTR_PXI_MEM_TYPE_BAR3           =  0x3FFF0214
const VI_ATTR_PXI_MEM_TYPE_BAR4           =  0x3FFF0215
const VI_ATTR_PXI_MEM_TYPE_BAR5           =  0x3FFF0216
const VI_ATTR_PXI_MEM_BASE_BAR0_32        =  0x3FFF0221
const VI_ATTR_PXI_MEM_BASE_BAR1_32        =  0x3FFF0222
const VI_ATTR_PXI_MEM_BASE_BAR2_32        =  0x3FFF0223
const VI_ATTR_PXI_MEM_BASE_BAR3_32        =  0x3FFF0224
const VI_ATTR_PXI_MEM_BASE_BAR4_32        =  0x3FFF0225
const VI_ATTR_PXI_MEM_BASE_BAR5_32        =  0x3FFF0226
const VI_ATTR_PXI_MEM_BASE_BAR0_64        =  0x3FFF0228
const VI_ATTR_PXI_MEM_BASE_BAR1_64        =  0x3FFF0229
const VI_ATTR_PXI_MEM_BASE_BAR2_64        =  0x3FFF022A
const VI_ATTR_PXI_MEM_BASE_BAR3_64        =  0x3FFF022B
const VI_ATTR_PXI_MEM_BASE_BAR4_64        =  0x3FFF022C
const VI_ATTR_PXI_MEM_BASE_BAR5_64        =  0x3FFF022D
const VI_ATTR_PXI_MEM_SIZE_BAR0_32        =  0x3FFF0231
const VI_ATTR_PXI_MEM_SIZE_BAR1_32        =  0x3FFF0232
const VI_ATTR_PXI_MEM_SIZE_BAR2_32        =  0x3FFF0233
const VI_ATTR_PXI_MEM_SIZE_BAR3_32        =  0x3FFF0234
const VI_ATTR_PXI_MEM_SIZE_BAR4_32        =  0x3FFF0235
const VI_ATTR_PXI_MEM_SIZE_BAR5_32        =  0x3FFF0236
const VI_ATTR_PXI_MEM_SIZE_BAR0_64        =  0x3FFF0238
const VI_ATTR_PXI_MEM_SIZE_BAR1_64        =  0x3FFF0239
const VI_ATTR_PXI_MEM_SIZE_BAR2_64        =  0x3FFF023A
const VI_ATTR_PXI_MEM_SIZE_BAR3_64        =  0x3FFF023B
const VI_ATTR_PXI_MEM_SIZE_BAR4_64        =  0x3FFF023C
const VI_ATTR_PXI_MEM_SIZE_BAR5_64        =  0x3FFF023D
const VI_ATTR_PXI_IS_EXPRESS              =  0x3FFF0240
const VI_ATTR_PXI_SLOT_LWIDTH             =  0x3FFF0241
const VI_ATTR_PXI_MAX_LWIDTH              =  0x3FFF0242
const VI_ATTR_PXI_ACTUAL_LWIDTH           =  0x3FFF0243
const VI_ATTR_PXI_DSTAR_BUS               =  0x3FFF0244
const VI_ATTR_PXI_DSTAR_SET               =  0x3FFF0245
const VI_ATTR_PXI_ALLOW_WRITE_COMBINE     =  0x3FFF0246
const VI_ATTR_TCPIP_HISLIP_OVERLAP_EN     =  0x3FFF0300
const VI_ATTR_TCPIP_HISLIP_VERSION        =  0x3FFF0301
const VI_ATTR_TCPIP_HISLIP_MAX_MESSAGE_KB =  0x3FFF0302
const VI_ATTR_TCPIP_IS_HISLIP             =  0x3FFF0303

const VI_ATTR_JOB_ID                      =  0x3FFF4006
const VI_ATTR_EVENT_TYPE                  =  0x3FFF4010
const VI_ATTR_SIGP_STATUS_ID              =  0x3FFF4011
const VI_ATTR_RECV_TRIG_ID                =  0x3FFF4012
const VI_ATTR_INTR_STATUS_ID              =  0x3FFF4023
const VI_ATTR_STATUS                      =  0x3FFF4025
const VI_ATTR_RET_COUNT_32                =  0x3FFF4026
const VI_ATTR_BUFFER                      =  0x3FFF4027
const VI_ATTR_RECV_INTR_LEVEL             =  0x3FFF4041
const VI_ATTR_OPER_NAME                   =  0xBFFF4042
const VI_ATTR_GPIB_RECV_CIC_STATE         =  0x3FFF4193
const VI_ATTR_RECV_TCPIP_ADDR             =  0xBFFF4198
const VI_ATTR_USB_RECV_INTR_SIZE          =  0x3FFF41B0
const VI_ATTR_USB_RECV_INTR_DATA          =  0xBFFF41B1
const VI_ATTR_PXI_RECV_INTR_SEQ           =  0x3FFF4240
const VI_ATTR_PXI_RECV_INTR_DATA          =  0x3FFF4241

#- Attributes (platform dependent size) ------------------------------------*/

if WORD_SIZE == 64
  ViAttrState = ViUInt64
else
  ViAttrState = ViUInt32
end

if WORD_SIZE == 64
  const VI_ATTR_USER_DATA_64                =  0x3FFF000A
  const VI_ATTR_RET_COUNT_64                =  0x3FFF4028
  const VI_ATTR_USER_DATA                   = VI_ATTR_USER_DATA_64
  const VI_ATTR_RET_COUNT                   = VI_ATTR_RET_COUNT_64
else
  const VI_ATTR_USER_DATA                   = VI_ATTR_USER_DATA_32
  const VI_ATTR_RET_COUNT                   = VI_ATTR_RET_COUNT_32
end

if WORD_SIZE == 64
  const VI_ATTR_WIN_BASE_ADDR_64            =  0x3FFF009B
  const VI_ATTR_WIN_SIZE_64                 =  0x3FFF009C
  const VI_ATTR_MEM_BASE_64                 =  0x3FFF00D0
  const VI_ATTR_MEM_SIZE_64                 =  0x3FFF00D1
end

if WORD_SIZE == 64
  const VI_ATTR_WIN_BASE_ADDR         = VI_ATTR_WIN_BASE_ADDR_64
  const VI_ATTR_WIN_SIZE              = VI_ATTR_WIN_SIZE_64
  const VI_ATTR_MEM_BASE              = VI_ATTR_MEM_BASE_64
  const VI_ATTR_MEM_SIZE              = VI_ATTR_MEM_SIZE_64
  const VI_ATTR_PXI_MEM_BASE_BAR0     = VI_ATTR_PXI_MEM_BASE_BAR0_64
  const VI_ATTR_PXI_MEM_BASE_BAR1     = VI_ATTR_PXI_MEM_BASE_BAR1_64
  const VI_ATTR_PXI_MEM_BASE_BAR2     = VI_ATTR_PXI_MEM_BASE_BAR2_64
  const VI_ATTR_PXI_MEM_BASE_BAR3     = VI_ATTR_PXI_MEM_BASE_BAR3_64
  const VI_ATTR_PXI_MEM_BASE_BAR4     = VI_ATTR_PXI_MEM_BASE_BAR4_64
  const VI_ATTR_PXI_MEM_BASE_BAR5     = VI_ATTR_PXI_MEM_BASE_BAR5_64
  const VI_ATTR_PXI_MEM_SIZE_BAR0     = VI_ATTR_PXI_MEM_SIZE_BAR0_64
  const VI_ATTR_PXI_MEM_SIZE_BAR1     = VI_ATTR_PXI_MEM_SIZE_BAR1_64
  const VI_ATTR_PXI_MEM_SIZE_BAR2     = VI_ATTR_PXI_MEM_SIZE_BAR2_64
  const VI_ATTR_PXI_MEM_SIZE_BAR3     = VI_ATTR_PXI_MEM_SIZE_BAR3_64
  const VI_ATTR_PXI_MEM_SIZE_BAR4     = VI_ATTR_PXI_MEM_SIZE_BAR4_64
  const VI_ATTR_PXI_MEM_SIZE_BAR5     = VI_ATTR_PXI_MEM_SIZE_BAR5_64
else
  const VI_ATTR_WIN_BASE_ADDR         = VI_ATTR_WIN_BASE_ADDR_32
  const VI_ATTR_WIN_SIZE              = VI_ATTR_WIN_SIZE_32
  const VI_ATTR_MEM_BASE              = VI_ATTR_MEM_BASE_32
  const VI_ATTR_MEM_SIZE              = VI_ATTR_MEM_SIZE_32
  const VI_ATTR_PXI_MEM_BASE_BAR0     = VI_ATTR_PXI_MEM_BASE_BAR0_32
  const VI_ATTR_PXI_MEM_BASE_BAR1     = VI_ATTR_PXI_MEM_BASE_BAR1_32
  const VI_ATTR_PXI_MEM_BASE_BAR2     = VI_ATTR_PXI_MEM_BASE_BAR2_32
  const VI_ATTR_PXI_MEM_BASE_BAR3     = VI_ATTR_PXI_MEM_BASE_BAR3_32
  const VI_ATTR_PXI_MEM_BASE_BAR4     = VI_ATTR_PXI_MEM_BASE_BAR4_32
  const VI_ATTR_PXI_MEM_BASE_BAR5     = VI_ATTR_PXI_MEM_BASE_BAR5_32
  const VI_ATTR_PXI_MEM_SIZE_BAR0     = VI_ATTR_PXI_MEM_SIZE_BAR0_32
  const VI_ATTR_PXI_MEM_SIZE_BAR1     = VI_ATTR_PXI_MEM_SIZE_BAR1_32
  const VI_ATTR_PXI_MEM_SIZE_BAR2     = VI_ATTR_PXI_MEM_SIZE_BAR2_32
  const VI_ATTR_PXI_MEM_SIZE_BAR3     = VI_ATTR_PXI_MEM_SIZE_BAR3_32
  const VI_ATTR_PXI_MEM_SIZE_BAR4     = VI_ATTR_PXI_MEM_SIZE_BAR4_32
  const VI_ATTR_PXI_MEM_SIZE_BAR5     = VI_ATTR_PXI_MEM_SIZE_BAR5_32
end

#- Event Types -------------------------------------------------------------*/

const VI_EVENT_IO_COMPLETION              =  0x3FFF2009
const VI_EVENT_TRIG                       =  0xBFFF200A
const VI_EVENT_SERVICE_REQ                =  0x3FFF200B
const VI_EVENT_CLEAR                      =  0x3FFF200D
const VI_EVENT_EXCEPTION                  =  0xBFFF200E
const VI_EVENT_GPIB_CIC                   =  0x3FFF2012
const VI_EVENT_GPIB_TALK                  =  0x3FFF2013
const VI_EVENT_GPIB_LISTEN                =  0x3FFF2014
const VI_EVENT_VXI_VME_SYSFAIL            =  0x3FFF201D
const VI_EVENT_VXI_VME_SYSRESET           =  0x3FFF201E
const VI_EVENT_VXI_SIGP                   =  0x3FFF2020
const VI_EVENT_VXI_VME_INTR               =  0xBFFF2021
const VI_EVENT_PXI_INTR                   =  0x3FFF2022
const VI_EVENT_TCPIP_CONNECT              =  0x3FFF2036
const VI_EVENT_USB_INTR                   =  0x3FFF2037

const VI_ALL_ENABLED_EVENTS               =  0x3FFF7FFF

#- Other VISA Definitions --------------------------------------------------*/

# const VI_VERSION_MAJOR(ver)       ((((ViVersion)ver) & 0xFFF00000 >> 20)
# const VI_VERSION_MINOR(ver)       ((((ViVersion)ver) & 0x000FFF00 >>  8)
# const VI_VERSION_SUBMINOR(ver)    ((((ViVersion)ver) & 0x000000FF      )

const VI_FIND_BUFLEN          = 256

const VI_INTF_GPIB            = 1
const VI_INTF_VXI             = 2
const VI_INTF_GPIB_VXI        = 3
const VI_INTF_ASRL            = 4
const VI_INTF_PXI             = 5
const VI_INTF_TCPIP           = 6
const VI_INTF_USB             = 7

const VI_PROT_NORMAL          = 1
const VI_PROT_FDC             = 2
const VI_PROT_HS488           = 3
const VI_PROT_4882_STRS       = 4
const VI_PROT_USBTMC_VENDOR   = 5

const VI_FDC_NORMAL           = 1
const VI_FDC_STREAM           = 2

const VI_LOCAL_SPACE          = 0
const VI_A16_SPACE            = 1
const VI_A24_SPACE            = 2
const VI_A32_SPACE            = 3
const VI_A64_SPACE            = 4
const VI_PXI_ALLOC_SPACE      = 9
const VI_PXI_CFG_SPACE        = 10
const VI_PXI_BAR0_SPACE       = 11
const VI_PXI_BAR1_SPACE       = 12
const VI_PXI_BAR2_SPACE       = 13
const VI_PXI_BAR3_SPACE       = 14
const VI_PXI_BAR4_SPACE       = 15
const VI_PXI_BAR5_SPACE       = 16
const VI_OPAQUE_SPACE         = 0xFFFF

const VI_UNKNOWN_LA               = -1
const VI_UNKNOWN_SLOT             = -1
const VI_UNKNOWN_LEVEL            = -1
const VI_UNKNOWN_CHASSIS          = -1

const VI_QUEUE                    = 1
const VI_HNDLR                    = 2
const VI_SUSPEND_HNDLR            = 4
const VI_ALL_MECH               =  0xFFFF

const VI_ANY_HNDLR                = 0

const VI_TRIG_ALL                 = -2
const VI_TRIG_SW                  = -1
const VI_TRIG_TTL0                = 0
const VI_TRIG_TTL1                = 1
const VI_TRIG_TTL2                = 2
const VI_TRIG_TTL3                = 3
const VI_TRIG_TTL4                = 4
const VI_TRIG_TTL5                = 5
const VI_TRIG_TTL6                = 6
const VI_TRIG_TTL7                = 7
const VI_TRIG_ECL0                = 8
const VI_TRIG_ECL1                = 9
const VI_TRIG_ECL2                = 10
const VI_TRIG_ECL3                = 11
const VI_TRIG_ECL4                = 12
const VI_TRIG_ECL5                = 13
const VI_TRIG_STAR_SLOT1          = 14
const VI_TRIG_STAR_SLOT2          = 15
const VI_TRIG_STAR_SLOT3          = 16
const VI_TRIG_STAR_SLOT4          = 17
const VI_TRIG_STAR_SLOT5          = 18
const VI_TRIG_STAR_SLOT6          = 19
const VI_TRIG_STAR_SLOT7          = 20
const VI_TRIG_STAR_SLOT8          = 21
const VI_TRIG_STAR_SLOT9          = 22
const VI_TRIG_STAR_SLOT10         = 23
const VI_TRIG_STAR_SLOT11         = 24
const VI_TRIG_STAR_SLOT12         = 25
const VI_TRIG_STAR_INSTR          = 26
const VI_TRIG_PANEL_IN            = 27
const VI_TRIG_PANEL_OUT           = 28
const VI_TRIG_STAR_VXI0           = 29
const VI_TRIG_STAR_VXI1           = 30
const VI_TRIG_STAR_VXI2           = 31

const VI_TRIG_PROT_DEFAULT        = 0
const VI_TRIG_PROT_ON             = 1
const VI_TRIG_PROT_OFF            = 2
const VI_TRIG_PROT_SYNC           = 5
const VI_TRIG_PROT_RESERVE        = 6
const VI_TRIG_PROT_UNRESERVE      = 7

const VI_READ_BUF                 = 1
const VI_WRITE_BUF                = 2
const VI_READ_BUF_DISCARD         = 4
const VI_WRITE_BUF_DISCARD        = 8
const VI_IO_IN_BUF                = 16
const VI_IO_OUT_BUF               = 32
const VI_IO_IN_BUF_DISCARD        = 64
const VI_IO_OUT_BUF_DISCARD       = 128

const VI_FLUSH_ON_ACCESS          = 1
const VI_FLUSH_WHEN_FULL          = 2
const VI_FLUSH_DISABLE            = 3

const VI_NMAPPED                  = 1
const VI_USE_OPERS                = 2
const VI_DEREF_ADDR               = 3
const VI_DEREF_ADDR_BYTE_SWAP     = 4

const VI_TMO_IMMEDIATE            = 0x00000000
const VI_TMO_INFINITE             = 0xFFFFFFFF

const VI_NO_LOCK                  = 0x00000000
const VI_EXCLUSIVE_LOCK           = 0x00000001
const VI_SHARED_LOCK              = 0x00000002
const VI_LOAD_CONFIG              = 0x00000003

const VI_NO_SEC_ADDR              = 0xFFFF

const VI_ASRL_PAR_NONE            = 0
const VI_ASRL_PAR_ODD             = 1
const VI_ASRL_PAR_EVEN            = 2
const VI_ASRL_PAR_MARK            = 3
const VI_ASRL_PAR_SPACE           = 4

const VI_ASRL_STOP_ONE            = 10
const VI_ASRL_STOP_ONE5           = 15
const VI_ASRL_STOP_TWO            = 20

const VI_ASRL_FLOW_NONE           = 0
const VI_ASRL_FLOW_XON_XOFF       = 1
const VI_ASRL_FLOW_RTS_CTS        = 2
const VI_ASRL_FLOW_DTR_DSR        = 4

const VI_ASRL_END_NONE            = 0
const VI_ASRL_END_LAST_BIT        = 1
const VI_ASRL_END_TERMCHAR        = 2
const VI_ASRL_END_BREAK           = 3

const VI_STATE_ASSERTED           = 1
const VI_STATE_UNASSERTED         = 0
const VI_STATE_UNKNOWN            = -1

const VI_BIG_ENDIAN               = 0
const VI_LITTLE_ENDIAN            = 1

const VI_DATA_PRIV                = 0
const VI_DATA_NPRIV               = 1
const VI_PROG_PRIV                = 2
const VI_PROG_NPRIV               = 3
const VI_BLCK_PRIV                = 4
const VI_BLCK_NPRIV               = 5
const VI_D64_PRIV                 = 6
const VI_D64_NPRIV                = 7
const VI_D64_2EVME                = 8
const VI_D64_SST160               = 9
const VI_D64_SST267               = 10
const VI_D64_SST320               = 11

const VI_WIDTH_8                  = 1
const VI_WIDTH_16                 = 2
const VI_WIDTH_32                 = 4
const VI_WIDTH_64                 = 8

const VI_GPIB_REN_DEASSERT        = 0
const VI_GPIB_REN_ASSERT          = 1
const VI_GPIB_REN_DEASSERT_GTL    = 2
const VI_GPIB_REN_ASSERT_ADDRESS  = 3
const VI_GPIB_REN_ASSERT_LLO      = 4
const VI_GPIB_REN_ASSERT_ADDRESS_LLO = 5
const VI_GPIB_REN_ADDRESS_GTL     = 6

const VI_GPIB_ATN_DEASSERT        = 0
const VI_GPIB_ATN_ASSERT          = 1
const VI_GPIB_ATN_DEASSERT_HANDSHAKE = 2
const VI_GPIB_ATN_ASSERT_IMMEDIATE = 3

const VI_GPIB_HS488_DISABLED      = 0
const VI_GPIB_HS488_NIMPL         = -1

const VI_GPIB_UNADDRESSED         = 0
const VI_GPIB_TALKER              = 1
const VI_GPIB_LISTENER            = 2

const VI_VXI_CMD16              =  0x0200
const VI_VXI_CMD16_RESP16       =  0x0202
const VI_VXI_RESP16             =  0x0002
const VI_VXI_CMD32              =  0x0400
const VI_VXI_CMD32_RESP16       =  0x0402
const VI_VXI_CMD32_RESP32       =  0x0404
const VI_VXI_RESP32             =  0x0004

const VI_ASSERT_SIGNAL            = -1
const VI_ASSERT_USE_ASSIGNED      = 0
const VI_ASSERT_IRQ1              = 1
const VI_ASSERT_IRQ2              = 2
const VI_ASSERT_IRQ3              = 3
const VI_ASSERT_IRQ4              = 4
const VI_ASSERT_IRQ5              = 5
const VI_ASSERT_IRQ6              = 6
const VI_ASSERT_IRQ7              = 7

const VI_UTIL_ASSERT_SYSRESET     = 1
const VI_UTIL_ASSERT_SYSFAIL      = 2
const VI_UTIL_DEASSERT_SYSFAIL    = 3

const VI_VXI_CLASS_MEMORY         = 0
const VI_VXI_CLASS_EXTENDED       = 1
const VI_VXI_CLASS_MESSAGE        = 2
const VI_VXI_CLASS_REGISTER       = 3
const VI_VXI_CLASS_OTHER          = 4

const VI_PXI_ADDR_NONE            = 0
const VI_PXI_ADDR_MEM             = 1
const VI_PXI_ADDR_IO              = 2
const VI_PXI_ADDR_CFG             = 3

const VI_TRIG_UNKNOWN             = -1

const VI_PXI_LBUS_UNKNOWN         = -1
const VI_PXI_LBUS_NONE            = 0
const VI_PXI_LBUS_STAR_TRIG_BUS_0 = 1000
const VI_PXI_LBUS_STAR_TRIG_BUS_1 = 1001
const VI_PXI_LBUS_STAR_TRIG_BUS_2 = 1002
const VI_PXI_LBUS_STAR_TRIG_BUS_3 = 1003
const VI_PXI_LBUS_STAR_TRIG_BUS_4 = 1004
const VI_PXI_LBUS_STAR_TRIG_BUS_5 = 1005
const VI_PXI_LBUS_STAR_TRIG_BUS_6 = 1006
const VI_PXI_LBUS_STAR_TRIG_BUS_7 = 1007
const VI_PXI_LBUS_STAR_TRIG_BUS_8 = 1008
const VI_PXI_LBUS_STAR_TRIG_BUS_9 = 1009
const VI_PXI_STAR_TRIG_CONTROLLER = 1413


text:
struct StiffWrapper{T,S} <: AbstractMatrix{T}
    parent::S
    max::T
    min::T
    penalty::T
    translation::T
    @inline function StiffWrapper{T,S}(A::S,penalty) where {T,S}
        maxV=typemin(T)
        minV=typemax(T)
        for i in eachindex(A)
            v=A[i]
            isinf(v) && continue
            maxV=max(maxV,v)
            minV=min(minV,v)
        end
        if isinf(maxV) || isinf(minV)
            maxV=minV=zero(T)
        end
        trasl=-minV
        if iszero(maxV+trasl)
            trasl+=one(T) #regularize illposed stiff problems
        end
        new(A,maxV,minV,penalty,trasl)
    end
end
@inline StiffWrapper(A,penalty)=StiffWrapper{eltype(A),typeof(A)}(A,penalty)

import Base.size, Base.getindex, Base.IndexStyle

@inline IndexStyle(::StiffWrapper) = IndexCartesian()
@inline function size_parent(A::StiffWrapper)
    size(A.parent)
end
@inline function size(A::StiffWrapper)
    n,m=size_parent(A)
    (n+m,n+m)
end
@inline function getindex(A::StiffWrapper{T},i::Int,j::Int) where {T}
    @boundscheck checkbounds(A,i,j)
    @inbounds begin
    n,m=size_parent(A)
    if i<=n && j<=m
        return A.parent[i,j]+A.translation
    end
    if i==j+n || j==i+m
        return (A.max+A.translation)*A.penalty
    end
    if i>n && j>m
        v=A.parent[j-m,i-n]+A.translation
        if !isinf(v)
            return A.min+A.translation
        end
    end
    return typemax(T)
    end
end

ext: jl

max_stars_repo_name: francescoalemanno/SolveLAP


text:
# These are functions defined in base/math.jl.  They apply elementwise
# to array arguments.  This is implemented in operators.jl using the
# @vectorize_1arg macro.  The first element of the tuples below gives
# the function name.  The second element is an expression for dy/dx
# where y=f(x).  The gradient expressions should be valid for array
# inputs.  The last part gives the range for generating test cases.

math1arg = [
(:cbrt, :(1./(3.*abs2(y))), (-Inf,Inf)),
(:deg2rad, :(pi/180), (-Inf,Inf)),
(:exp, :y, (-Inf,Inf)),
(:exp10, :(y.*log(10)), (-Inf,Inf)),
(:exp2, :(y.*log(2)), (-Inf,Inf)),
(:expm1, :(1+y), (-Inf,Inf)),
(:log, :(1./x), (0,Inf)),
(:log10, :(1./(log(10).*x)), (0,Inf)),
(:log1p, :(1./(1+x)), (-1,Inf)),
(:log2, :(1./(log(2).*x)), (0,Inf)),
(:rad2deg, :(180/pi), (-Inf,Inf)),
(:significand, :(0.5.^exponent(x)), (-Inf,Inf)),
(:sqrt, :(1./(2.*y)), (0,Inf)),
]

for (f,g,r) in math1arg
    @eval @primitive $f(x),dy,y  (dy.*($g))
    addtest1(f,r)
end

# math2arg: These are two argument functions that can handle mixing
# scalar and array arguments.  They use vectorize_2arg in Base
# (defined in operators) which allows them to have Array arguments in
# first, second, or both positions.  When both arguments are Arrays
# they must have the same size, or if one has extra dimensions at the
# end, they need to be 1.  The resulting array will have the longer of
# the two sizes.  (implemented by promote_shape).  Note that no
# broadcasting is performed here, i.e. the two arrays need to have the
# same length.  Note that some 2arg functions are defined in arraymath
# not using vectorize_2arg.  Using variable names: y=f(x1,x2) in
# gradient definitions.

math2arg = [
(:atan2, :(x2./(abs2(x1)+abs2(x2))), :(-x1./(abs2(x1)+abs2(x2)))),
(:hypot, :(x1./y), :(x2./y)),
(:max, :(y.==x1), :(y.==x2)),
(:min, :(y.==x1), :(y.==x2)),
]

for (f,g1,g2) in math2arg
    @eval @primitive $f(x1,x2),dy,y  unbroadcast(x1,dy.*($g1))  unbroadcast(x2,dy.*($g2))
    addtest2(f,(-Inf,Inf))
end

# The 2-arg log supports positive args for reals.
log(x1::Irrational{:e},x2::Value)=log(float(x1),x2) # to avoid clash with irrationals.jl:131.
@primitive log(x1,x2),dy  unbroadcast(x1,-dy.*log(x2)./(x1.*abs2(log(x1))))  unbroadcast(x2,dy./(x2.*log(x1)))
addtest2(log,(0,Inf))

# ^ only supports (N>=0,N), arrays not supported in math.jl, only M^N in linalg/dense.jl (TODO)
(^){T<:Number}(x1::Value{T},x2::Integer)=(^)(x1,float(x2)) # to avoid clash with intfuncs:108
@primitive (^)(x1::Number,x2::Number),dy,y  (dy*x2*x1^(x2-1))  (dy*y*log(x1))
addtest(^, randin((0,Inf)), randin((-Inf,Inf)))

# clamp(x,lo,hi) clamps x between lo and hi
@primitive clamp(x,i...),dy,y  unbroadcast(x,dy.*(i[1] .<= x .<= i[2]))
addtest(clamp, randn(10), -1., 1.)
addtest(clamp, randn(), -1., 1.)

# ldexp(x,n) computes x*2^n with x real, n integer
@primitive ldexp(x,n...),dy  (dy*(2.0^n[1]))
addtest(ldexp, randn(), rand(-2:2))

# mod2pi(x) returns modulus after division by 2pi for x real.
@primitive mod2pi(x::Number),dy dy
addtest(mod2pi, 100randn())

# zerograd functions
@zerograd exponent(x)

# Other functions defined in julia/base/math.jl
# add22condh: Not exported
# angle_restrict_symm: Not exported
# clamp!: overwriting
# frexp: returns tuple
# ieee754_rem_pio2: Not exported
# minmax: returns tuple
# modf: returns tuple
# Moved to erf.jl:
# erf: see erf.jl
# erfc: see erf.jl
# The following moved to trig.jl:
# (acos, :(-1./sqrt(1-abs2(x))), (-1,1)),
# (acosh, :(1./sqrt(abs2(x)-1)), (1,Inf)),
# (asin, :(1./sqrt(1-abs2(x))), (-1,1)),
# (asinh, :(1./sqrt(1+abs2(x))), (-Inf,Inf)),
# (atan, :(1./(1+abs2(x))), (-Inf,Inf)),
# (atanh, :(1./(1-abs2(x))), (-1,1)),
# (cos, :(-sin(x)), (-Inf,Inf)),
# (cosh, :(sinh(x)), (-Inf,Inf)),
# (sin, :(cos(x)), (-Inf,Inf)),
# (sinh, :(cosh(x)), (-Inf,Inf)),
# (tan, :(1+abs2(y)), (-Inf,Inf)),
# (tanh, :(1-abs2(y)), (-Inf,Inf)),
# Moved to gamma.jl:
# (lgamma, :(digamma(x)), (-Inf,Inf)),


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 728d5546-7a96-11eb-1757-93681c4e7955
begin
	import Pkg
    Pkg.activate(".")
    Pkg.add("PlutoUI")
	Pkg.add("Plots")

	using PlutoUI
	using Plots
end

# ╔═╡ 58ba05d0-7a96-11eb-08da-41b782596598
md"# Euler Method"

# ╔═╡ 97f8a9f8-7b78-11eb-2f25-e5228f41db56
md"We implement a Euer Method to solve first oderd ODEs and use it to some relevant Biological Models."

# ╔═╡ ecad4e78-7a96-11eb-2b54-e7f910e83835
begin
	function euler(f,t0,tn,Δt,x₀)
		T = t0:Δt:tn
		F = Vector{typeof(x₀)}(undef,length(T))
		F[1] = x₀
		for (n,t) ∈ enumerate(T[2:end])
			F[n+1] = F[n] .+ f(F[n]) .* Δt
		end
		return F
	end

	close(f,par) = x -> f(x,par...)

	reshape_result(F) = [[f[i] for f ∈ F] for i ∈ 1:length(F[1])] ;
end

# ╔═╡ cd0f0650-7b78-11eb-12c0-cb427957850b
md"""
---
### Logistic Population Growth
"""

# ╔═╡ e9b52262-7b7d-11eb-21aa-2bf0d9f5ef50
md"""

$\frac{dP}{dt} = rP \left( 1- \frac{P}{K} \right)$

"""

# ╔═╡ dcffb046-7b78-11eb-333c-7df614c12617
begin
	LogisticGrowth(x,K,r) = r*x*(1-x/K)
	LG_0 = 1.0
end;

# ╔═╡ 04fdd000-7b79-11eb-029b-b3765364f627
md"""
growth rate r $(@bind r NumberField(0.0001:0.0001:0.1,default=0.01))

carrying capacity K $(@bind K NumberField(1:10^7,default=10_000))
"""

# ╔═╡ 5b3c4954-7b79-11eb-11ba-57b8ae51a6f5
let
	Δt = 0.01

	p = plot(framestlye=:zerolines,legend=:topleft)

	T=0:Δt:2500
	F = euler(close(LogisticGrowth,[K,r]),T[1],T[end],Δt,LG_0)
	plot!(p,T,F,label="Population Size")

	p
end

# ╔═╡ 7a010972-7b78-11eb-3478-f1fbc1f0645d
md"""
---
### Lotka-Volterra
"""

# ╔═╡ b963fdb2-7b7e-11eb-3184-1b216f77bd64
md"""

${\frac{dN_1}{dt} = N_1 (\varepsilon_1 - \gamma_1 N_2) \quad , \quad \frac{dN_2}{dt} = -N_2 (\varepsilon_2 - \gamma_2 N_1)}$

"""

# ╔═╡ 7fc751be-7a96-11eb-207d-d5d0142f43ba
begin
	LotkaVolterra(x,ε₁,δ₁,ε₂,δ₂) = [
			x[1]*(ε₁-δ₁*x[2]),
			-x[2]*(ε₂-δ₂*x[1])
	]
	LV_0 = [10.0,10.0]
end;

# ╔═╡ 07b392ee-7b77-11eb-00be-876772164be8
md"""
Reproduktionsraten
	ε₁ $(@bind ε₁ Slider(0.01:0.01:1,show_value=true, default=0.5))
	ε₂ $(@bind ε₂ Slider(0.01:0.01:1,show_value=true, default=0.25))

Sterberaten
	δ₁ $(@bind δ₁ Slider(0.01:0.01:1,show_value=true, default=0.07))
	δ₂ $(@bind δ₂ Slider(0.01:0.01:1,show_value=true, default=0.07))
"""

# ╔═╡ c080dd1a-7a96-11eb-345c-e364dc0f8ac9
let
	Δt = 0.01

	p = plot(framestlye=:zerolines)

	T=0:Δt:100
	F = euler(close(LotkaVolterra,[ε₁,δ₁,ε₂,δ₂]),T[1],T[end],Δt,LV_0)
	plot!(p,T,reshape_result(F),label=["N₁" "N₂"])

	p
end

# ╔═╡ 62645b8c-7b7a-11eb-1b35-a3957a0252fb
md"""
---
### SIR-Model
"""

# ╔═╡ 4341adf4-7b7f-11eb-25c1-e546febcad3e
md"""

${\frac{dS}{dt} = \nu N - \beta \frac{SI}{N} - \mu S}$

${\frac{dI}{dt} = \beta \frac{SI}{N} - \gamma I - \mu I}$

${\frac{dR}{dt} = \gamma I - \mu R \phantom{-\beta \frac{SI}{N}}}$

"""

# ╔═╡ 6b37eb20-7b7a-11eb-07de-a1042efa6721
begin
	function SIR(x,ν,β,μ,γ)
		N = sum(x)
		S, I, R = x[1], x[2], x[3]
		return [
			(ν*N - β*S*I/N - μ*S),	# dS/dt
			(β*S*I/N - γ*I - μ*I),	# dI/dt
			(γ*I - μ*R)				# dR/dt
		]
	end

	SIR_0 = [9999.9,0.1,0.0]

	Δt = 0.01
	T=0:Δt:150
end;

# ╔═╡ e64e9d62-7b7b-11eb-3620-a7a57b954fcc
md"""
recovery rate γ $(@bind γ Slider(0.0:0.01:1,show_value=true,default=0.15)) |
infection rate β $(@bind β Slider(0.0:0.01:1,show_value=true,default=0.7))

death rate μ $(@bind μ Slider(0.0:0.01:1,show_value=true)) |
birth rate ν $(@bind ν Slider(0.0:0.01:1,show_value=true))
"""

# ╔═╡ fa94e89e-7c1e-11eb-08a7-91210bc2fe60
begin
	F = euler(close(SIR,[ν,β,μ,γ]),T[1],T[end],Δt,SIR_0)
	rF = reshape_result(F)
	N = sum.(F)
end;

# ╔═╡ 68f66db2-7b7c-11eb-138d-2f40b3bfbb4e
let
	p = plot(framestlye=:zerolines,size=(680,450))

	plot!(p,T,rF,label=["S" "I" "R"])

	p
end

# ╔═╡ Cell order:
# ╟─58ba05d0-7a96-11eb-08da-41b782596598
# ╟─97f8a9f8-7b78-11eb-2f25-e5228f41db56
# ╠═ecad4e78-7a96-11eb-2b54-e7f910e83835
# ╟─cd0f0650-7b78-11eb-12c0-cb427957850b
# ╟─e9b52262-7b7d-11eb-21aa-2bf0d9f5ef50
# ╠═dcffb046-7b78-11eb-333c-7df614c12617
# ╟─04fdd000-7b79-11eb-029b-b3765364f627
# ╠═5b3c4954-7b79-11eb-11ba-57b8ae51a6f5
# ╟─7a010972-7b78-11eb-3478-f1fbc1f0645d
# ╟─b963fdb2-7b7e-11eb-3184-1b216f77bd64
# ╠═7fc751be-7a96-11eb-207d-d5d0142f43ba
# ╟─07b392ee-7b77-11eb-00be-876772164be8
# ╠═c080dd1a-7a96-11eb-345c-e364dc0f8ac9
# ╟─62645b8c-7b7a-11eb-1b35-a3957a0252fb
# ╟─4341adf4-7b7f-11eb-25c1-e546febcad3e
# ╠═6b37eb20-7b7a-11eb-07de-a1042efa6721
# ╠═fa94e89e-7c1e-11eb-08a7-91210bc2fe60
# ╟─e64e9d62-7b7b-11eb-3620-a7a57b954fcc
# ╟─68f66db2-7b7c-11eb-138d-2f40b3bfbb4e
# ╟─728d5546-7a96-11eb-1757-93681c4e7955

max_stars_repo_name: roccminton/PlutoNotebooks

max_issues_repo_name: roccminton/PlutoNotebooks

max_forks_repo_name: roccminton/PlutoNotebooks


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# (This is part of the FFTW module.)

export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!

# Discrete cosine transforms (type II/III) via FFTW's r2r transforms;
# we follow the Matlab convention and adopt a unitary normalization here.
# Unlike Matlab we compute the multidimensional transform by default,
# similar to the Julia fft functions.

mutable struct DCTPlan{T<:fftwNumber,K,inplace} <: Plan{T}
    plan::r2rFFTWPlan{T}
    r::Array{UnitRange{Int}} # array of indices for rescaling
    nrm::Float64 # normalization factor
    region::Dims # dimensions being transformed
    pinv::DCTPlan{T}
    DCTPlan{T,K,inplace}(plan,r,nrm,region) where {T<:fftwNumber,K,inplace} = new(plan,r,nrm,region)
end

size(p::DCTPlan) = size(p.plan)

function show{T,K,inplace}(io::IO, p::DCTPlan{T,K,inplace})
    print(io, inplace ? "FFTW in-place " : "FFTW ",
          K == REDFT10 ? "DCT (DCT-II)" : "IDCT (DCT-III)", " plan for ")
    showfftdims(io, p.plan.sz, p.plan.istride, eltype(p))
end

for (pf, pfr, K, inplace) in ((:plan_dct, :plan_r2r, REDFT10, false),
                              (:plan_dct!, :plan_r2r!, REDFT10, true),
                              (:plan_idct, :plan_r2r, REDFT01, false),
                              (:plan_idct!, :plan_r2r!, REDFT01, true))
    @eval function $pf{T<:fftwNumber}(X::StridedArray{T}, region; kws...)
        r = [1:n for n in size(X)]
        nrm = sqrt(0.5^length(region) * normalization(X,region))
        DCTPlan{T,$K,$inplace}($pfr(X, $K, region; kws...), r, nrm,
                               ntuple(i -> Int(region[i]), length(region)))
    end
end

"""
    plan_dct!(A [, dims [, flags [, timelimit]]])

Same as [`plan_dct`](@ref), but operates in-place on `A`.
"""
plan_dct!

"""
    plan_idct(A [, dims [, flags [, timelimit]]])

Pre-plan an optimized inverse discrete cosine transform (DCT), similar to
[`plan_fft`](@ref) except producing a function that computes
[`idct`](@ref). The first two arguments have the same meaning as for
[`idct`](@ref).
"""
plan_idct

"""
    plan_dct(A [, dims [, flags [, timelimit]]])

Pre-plan an optimized discrete cosine transform (DCT), similar to
[`plan_fft`](@ref) except producing a function that computes
[`dct`](@ref). The first two arguments have the same meaning as for
[`dct`](@ref).
"""
plan_dct

"""
    plan_idct!(A [, dims [, flags [, timelimit]]])

Same as [`plan_idct`](@ref), but operates in-place on `A`.
"""
plan_idct!

function plan_inv{T,K,inplace}(p::DCTPlan{T,K,inplace})
    X = Array{T}(p.plan.sz)
    iK = inv_kind[K]
    DCTPlan{T,iK,inplace}(inplace ?
                          plan_r2r!(X, iK, p.region, flags=p.plan.flags) :
                          plan_r2r(X, iK, p.region, flags=p.plan.flags),
                          p.r, p.nrm, p.region)
end

for f in (:dct, :dct!, :idct, :idct!)
    pf = Symbol("plan_", f)
    @eval begin
        $f(x::AbstractArray{<:fftwNumber}) = $pf(x) * x
        $f(x::AbstractArray{<:fftwNumber}, region) = $pf(x, region) * x
        $pf(x::AbstractArray; kws...) = $pf(x, 1:ndims(x); kws...)
        $f(x::AbstractArray{<:Real}, region=1:ndims(x)) = $f(fftwfloat(x), region)
        $pf(x::AbstractArray{<:Real}, region; kws...) = $pf(fftwfloat(x), region; kws...)
        $pf(x::AbstractArray{<:Complex}, region; kws...) = $pf(fftwcomplex(x), region; kws...)
    end
end

"""
    dct(A [, dims])

Performs a multidimensional type-II discrete cosine transform (DCT) of the array `A`, using
the unitary normalization of the DCT. The optional `dims` argument specifies an iterable
subset of dimensions (e.g. an integer, range, tuple, or array) to transform along.  Most
efficient if the size of `A` along the transformed dimensions is a product of small primes;
see [`nextprod`](@ref). See also [`plan_dct`](@ref) for even greater
efficiency.
"""
dct

"""
    idct(A [, dims])

Computes the multidimensional inverse discrete cosine transform (DCT) of the array `A`
(technically, a type-III DCT with the unitary normalization). The optional `dims` argument
specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to
transform along.  Most efficient if the size of `A` along the transformed dimensions is a
product of small primes; see [`nextprod`](@ref).  See also
[`plan_idct`](@ref) for even greater efficiency.
"""
idct

"""
    dct!(A [, dims])

Same as [`dct!`](@ref), except that it operates in-place on `A`, which must be an
array of real or complex floating-point values.
"""
dct!

"""
    idct!(A [, dims])

Same as [`idct!`](@ref), but operates in-place on `A`.
"""
idct!

const sqrthalf = sqrt(0.5)
const sqrt2 = sqrt(2.0)
const onerange = 1:1

function A_mul_B!{T}(y::StridedArray{T}, p::DCTPlan{T,REDFT10},
                     x::StridedArray{T})
    assert_applicable(p.plan, x, y)
    unsafe_execute!(p.plan, x, y)
    scale!(y, p.nrm)
    r = p.r
    for d in p.region
        oldr = r[d]
        r[d] = onerange
        y[r...] *= sqrthalf
        r[d] = oldr
    end
    return y
end

# note: idct changes input data
function A_mul_B!{T}(y::StridedArray{T}, p::DCTPlan{T,REDFT01},
                     x::StridedArray{T})
    assert_applicable(p.plan, x, y)
    scale!(x, p.nrm)
    r = p.r
    for d in p.region
        oldr = r[d]
        r[d] = onerange
        x[r...] *= sqrt2
        r[d] = oldr
    end
    unsafe_execute!(p.plan, x, y)
    return y
end

*{T}(p::DCTPlan{T,REDFT10,false}, x::StridedArray{T}) =
    A_mul_B!(Array{T}(p.plan.osz), p, x)

*{T}(p::DCTPlan{T,REDFT01,false}, x::StridedArray{T}) =
    A_mul_B!(Array{T}(p.plan.osz), p, copy(x)) # need copy to preserve input

*{T,K}(p::DCTPlan{T,K,true}, x::StridedArray{T}) = A_mul_B!(x, p, x)

max_stars_repo_name: felixrehren/julia

max_issues_repo_name: felixrehren/julia

max_forks_repo_name: felixrehren/julia


text:
type CallSpecifiedMultiProduct{M <: MarketModelMultiProduct, E <: ExerciseStrategy, M2 <: MarketModelMultiProduct} <: MarketModelMultiProduct
  underlying::M
  strategy::E
  rebate::M2
  evolution::EvolutionDescription
  isPresent::Vector{BitArray{1}}
  cashFlowTimes::Vector{Float64}
  rebateOffset::Int
  wasCalled::Bool
  dummyCashFlowsThisStep::Vector{Int}
  dummyCashFlowsGenerated::Vector{Vector{MarketModelCashFlow}}
  currentIndex::Int
  callable::Bool
end

function CallSpecifiedMultiProduct(underlying::MarketModelMultiProduct,
                                  strategy::ExerciseStrategy,
                                  rebate::MarketModelMultiProduct)

  underlying = clone(underlying)
  strategy = clone(strategy)
  rebate = clone(rebate)

  products = number_of_products(underlying)
  d1 = get_evolution(underlying)
  rateTimes1 = d1.rateTimes
  evolutionTimes1 = d1.evolutionTimes
  exerciseTimes = strategy.exerciseTimes

  # TODO check if rebate isn't there
  d2 = get_evolution(rebate)
  rateTimes2 = d2.rateTimes
  rateTimes1 == rateTimes2 || error("incompatible rate times")

  allEvolutionTimes = Vector{Vector{Float64}}(4)
  allEvolutionTimes[1] = evolutionTimes1
  allEvolutionTimes[2] = exerciseTimes
  allEvolutionTimes[3] = d2.evolutionTimes
  allEvolutionTimes[4] = relevant_times(strategy)

  mergedEvolutionTimes, isPresent = merge_times(allEvolutionTimes)

  evolution = EvolutionDescription(rateTimes1, mergedEvolutionTimes)
  cashFlowTimes = possible_cash_flow_times(underlying)
  rebateOffset = length(cashFlowTimes)
  rebateTimes = possible_cash_flow_times(rebate)
  cashFlowTimes = vcat(cashFlowTimes, rebateTimes)

  dummyCashFlowsThisStep = zeros(Int, products)
  n = max_number_of_cashflows_per_product_per_step(rebate)

  dummyCashFlowsGenerated = fill(Vector{MarketModelCashFlow}(n), products)

  return CallSpecifiedMultiProduct(underlying, strategy, rebate, evolution, isPresent, cashFlowTimes, rebateOffset, false, dummyCashFlowsThisStep,
                                  dummyCashFlowsGenerated, 1, true)
end

get_evolution(cs::CallSpecifiedMultiProduct) = cs.evolution
number_of_products(cs::CallSpecifiedMultiProduct) = number_of_products(cs.underlying)
possible_cash_flow_times(cs::CallSpecifiedMultiProduct) = cs.cashFlowTimes
function reset!(cs::CallSpecifiedMultiProduct)
  reset!(cs.underlying)
  reset!(cs.rebate)
  reset!(cs.strategy)
  cs.currentIndex = 1
  cs.wasCalled = false

  return cs
end

max_number_of_cashflows_per_product_per_step(cs::CallSpecifiedMultiProduct) =
              max(max_number_of_cashflows_per_product_per_step(cs.underlying), max_number_of_cashflows_per_product_per_step(cs.rebate))

function next_time_step!(cs::CallSpecifiedMultiProduct, currentState::CurveState, numberCashFlowsThisStep::Vector{Int}, genCashFlows::Vector{Vector{MarketModelCashFlow}})
  isUnderlyingTime = cs.isPresent[1][cs.currentIndex]
  isExerciseTime = cs.isPresent[2][cs.currentIndex]
  isRebateTime = cs.isPresent[3][cs.currentIndex]
  isStrategyRelevantTime = cs.isPresent[4][cs.currentIndex]

  isDone = false

  if ~cs.wasCalled && isStrategyRelevantTime
    next_step!(cs.strategy, currentState)
  end

  if ~cs.wasCalled && isExerciseTime && cs.callable
    cs.wasCalled = get_exercise!(cs.strategy, currentState)
  end

  if cs.wasCalled
    if isRebateTime
      isDone = next_time_step!(cs.rebate, currentState, numberCashFlowsThisStep, genCashFlows)
      for i in eachindex(numberCashFlowsThisStep)
        for j = 1:numberCashFlowsThisStep[i]
          genCashFlows[i][j].timeIndex += cs.rebateOffset
        end
      end
    end
  else
    if isRebateTime
      next_time_step!(cs.rebate, currentState, cs.dummyCashFlowsThisStep, cs.dummyCashFlowsGenerated)
    end

    if isUnderlyingTime
      isDone = next_time_step!(cs.underlying, currentState, numberCashFlowsThisStep, genCashFlows)
    end
  end

  cs.currentIndex += 1
  return isDone || cs.currentIndex == length(cs.evolution.evolutionTimes) + 1
end

clone(cs::CallSpecifiedMultiProduct) = CallSpecifiedMultiProduct(clone(cs.underlying), clone(cs.strategy), clone(cs.rebate), clone(cs.evolution), deepcopy(cs.isPresent),
                                        copy(cs.cashFlowTimes), cs.rebateOffset, cs.wasCalled, copy(cs.dummyCashFlowsThisStep), deepcopy(cs.dummyCashFlowsGenerated),
                                        cs.currentIndex, cs.callable)


text:
################################################################################
function anchor_point(t,y,detail_th,split_wfm::Bool)
    if split_wfm
        nr_split = max(Int(round(length(y)/5e3)),1)
        #remove extra points for reshape TO BE SOLVED!!!
        aux = rem(length(y),nr_split)
        dt  = t[2] - t[1]
        y   = y[1:end-aux]
        t   = t[1:end-aux]
        t  = reshape(t,Int(length(t)/nr_split),nr_split)
        y  = reshape(y,Int(length(y)/nr_split),nr_split)
        # Collect the results as they become available.
        t_new    = [0.0, 0.0]
        y_new    = [0.0, 0.0]
        idx2save = [0, 0]
        for inr = 1:nr_split
            t_aux, y_aux, idx2save_aux = anchor_point(t[:,inr],y[:,inr],detail_th)
            append!(t_new,t_aux)
            #t_new      = t_new[1:end-1]
            append!(y_new,y_aux)
            #y_new      = y_new[1:end-1]
            append!(idx2save,idx2save_aux)
            #idx2save   = idx2save[1:end-1]
            #append!(idx2save2,idx2save_aux + (inr - 1)*Int64(length(y1)))TO BE SOLVED!!!
        end
        t_new = t_new[3:end]
        y_new = y_new[3:end]
    else
        t_new, y_new, idx2save = anchor_point(t,y,detail_th)
    end
    return t_new, y_new, idx2save
end
################################################################################
function anchor_point(t,y,detail_th)
    # Compress signal by removing redundant points.
    # Adjust waveform detail/compression ratio with 'detail_th' (maximum allowed
    # difference between original and approximated points from the signal)
    # For best performance perform wavelet noise reduction first
    yln           = length(y)
    idx_l         = 1
    idx_r         = yln
    idx_d_max     = 1
    cond_break    = true
    d_max         = 0.0
    M             = zeros(Int,yln) # hash table for relevant indices
    idx2save      = zeros(Int,yln+2)
    cnt           = 2
    idx2save[1:2] = [1,yln]
    while cond_break
        # get maximum error(difference) and index, between original chunk of signal
        # and linear approximation
        d_max, idx_d_max = get_d_max(idx_l,idx_r,y)
        # save all indices
        M[idx_d_max] = idx_r
        if d_max > detail_th
            # if computed error is greater than maximum allowed error, save
            # next point index and call get_d_max(idx_l,idx_r,y) at next
            # iteration; keep going towards leftmost branches
            cnt           = cnt + 1
            idx_r         = idx_d_max
            idx2save[cnt] = idx_d_max
        else
            # if computed error is smaller than maximum allowed error, stop, go
            # right(to the next waveform segment) and call get_d_max(idx_l,idx_r,y)
            # at the next iteration
            idx_l     = idx_r;
            if idx_l != yln
                idx_r = M[idx_l]
            else
                cond_break = false
            end
        end
    end
    # sort all indexes corresponding to relevent points and generate resampled
    # signal
    idx2save = idx2save[1:cnt]
    idx2save = sort(idx2save)
    t_new    = @view t[idx2save]
    y_new    = @view y[idx2save]
    return t_new, y_new, idx2save
end
################################################################################
function get_d_max(idx_l,idx_r,y)
    # cut segment to be resampled
    yp = view(y,idx_l:idx_r)
    # construct linear approximation
    dr = LinRange(y[idx_l], y[idx_r], length(yp))
    # compute distance(error) and get index of maximum error
    # -> this will be used for further splitting the
    # signal and will be part of the final resampled signal
    d_max     = 0.0
    idx_d_max = 1
    err_val   = 0.0
    for i = 1:length(yp)
        err_val = abs(yp[i] - dr[i])
        if err_val > d_max
            d_max     = err_val
            idx_d_max = i
        end
    end
    idx_d_max = idx_d_max + idx_l - 1
    return d_max, idx_d_max
end
################################################################################
# function anchor_point_gif(t,y,detail_th)
#     # Compress signal by removing redundant points.
#     # Adjust waveform detail/compression ratio with 'detail_th' (maximum allowed
#     # difference between original and approximated points from the signal)
#     # For best performance perform wavelet noise reduction first
#     # This version is written just for fun demonstration of the algorithm
#     yln           = length(y)
#     idx_l         = 1
#     idx_r         = yln
#     idx_d_max     = 1
#     cond_break    = true
#     d_max         = 0.0
#     M             = zeros(Int,yln) # hash table for relevant indices
#     idx2save      = zeros(Int,yln+2)
#     cnt           = 2
#     idx2save[1:2] = [1,yln]
#     anim = @animate for i = 1:(2*length(y))
#         # get maximum error(difference) and index, between original chunk of signal
#         # and linear approximation
#         d_max, idx_d_max = get_d_max(idx_l,idx_r,y)
#         # save all indices
#         M[idx_d_max] = idx_r
#         if d_max > detail_th
#             # if computed error is greater than maximum allowed error, save
#             # next point index and call get_d_max(idx_l,idx_r,y) at next
#             # iteration; keep going towards leftmost branches
#             cnt           = cnt + 1
#             idx_r         = idx_d_max
#             idx2save[cnt] = idx_d_max
#             idx2save_aux = sort(idx2save)
#             plt = plot(t,y,color=:blue, markershape=:circle, markersize=3, linealpha=0.2, legend=:none)
#             plot!(plt, t[idx2save_aux[idx2save_aux.!=0]],y[idx2save_aux[idx2save_aux.!=0]],
#             color=:red, markershape=:circle, markersize=5, linealpha=1, legend=:none)
#         else
#             # if computed error is smaller than maximum allowed error, stop, go
#             # right(to the next waveform segment) and call get_d_max(idx_l,idx_r,y)
#             # at the next iteration
#             idx_l     = idx_r;
#             if idx_l != yln
#                 idx_r = M[idx_l]
#             else
#                 break
#             end
#         end
#     end
#     # sort all indexes corresponding to relevent points and generate resampled
#     # signal
#     idx2save = idx2save[1:cnt]
#     idx2save = sort(idx2save)
#     t_new    = @view t[idx2save]
#     y_new    = @view y[idx2save]
#     return t_new, y_new, idx2save, anim
# end


text:
module TestPartiallyTypedTables

import Tables
using DataFrames: DataFrame, DataFrames
using PartiallyTypedTables: partiallytyped
using Test

# Monkey patch!
Tables.rows(df::DataFrames.AbstractDataFrame) = eachrow(df)

@testset begin
    df = DataFrame(a = [1, 3], b = [2.0, 4.0], c = [:a, :b], d = [nothing, nothing])
    @test map(row -> (row.a, row.b), partiallytyped(df, :a)) == [(1, 2), (3, 4)]

    demo(itr) = (first(itr).a, first(itr).b, first(itr).c, first(itr).d)
    @test @inferred(demo(partiallytyped(df, :a, :b, :c, :d))) == (1, 2.0, :a, nothing)
end

end  # module


text:
# This file contains functions for calling, running and post-processing
# SIIP PSI Simulation for actual Energy and A/S Market Clearing

"""
This function returns realized capacity factors for Thermal generators from PSI Simulation.
"""
function get_realized_capacity_factors(device::PSY.ThermalStandard,
                                        results::Dict{Symbol, DataFrames.DataFrame}
                                        )
    energy_production = results[:P__ThermalStandard][:, Symbol(get_name(device))]
    capacity_factors = energy_production / PSY.get_active_power_limits(device)[:max]
    return capacity_factors
end

"""
This function returns realized capacity factors for Renewable generators from PSI Simulation.
"""
function get_realized_capacity_factors(device::PSY.RenewableDispatch,
                                        results::Dict{Symbol, DataFrames.DataFrame}
                                        )
    energy_production = results[:P__RenewableDispatch][:, Symbol(get_name(device))]
    capacity_factors = energy_production / get_device_size(device)
    return capacity_factors
end

"""
This function returns realized capacity factors for Hydropower generators from PSI Simulation.
"""
function get_realized_capacity_factors(device::PSY.HydroEnergyReservoir,
                                        results::Dict{Symbol, DataFrames.DataFrame}
                                        )
    energy_production = results[:P__HydroEnergyReservoir][:, Symbol(get_name(device))]
    capacity_factors = energy_production / get_device_size(device)
    return capacity_factors
end

"""
This function returns nothing if Service is not of ReserveUp or ReserveDown type.
"""
function update_realized_reserve_perc!(device::PSY.Device,
                                        service::S,
                                        results::Dict{Symbol, DataFrames.DataFrame},
                                        reserve_up_perc::Dict{String, Array{Float64, 2}},
                                        reserve_down_perc::Dict{String, Array{Float64, 2}},
                                        data_length::Int64) where S <: PSY.Service
    return
end

"""
This function returns realized reserve up provision percentages from PSI Simulation.
"""
function update_realized_reserve_perc!(device::PSY.Device,
                                        service::PSY.VariableReserve{PSY.ReserveUp},
                                        results::Dict{Symbol, DataFrames.DataFrame},
                                        reserve_up_perc::Dict{String, Array{Float64, 2}},
                                        reserve_down_perc::Dict{String, Array{Float64, 2}})

    reserve_provision = results[Symbol("$(PSY.get_name(service))__VariableReserve_PowerSystems.ReserveUp")][:, Symbol(get_name(device))]
    reserve_perc = reserve_provision / get_device_size(device)
    reserve_up_perc[get_name(device)][1, :] = reserve_perc
    return
end

"""
This function returns realized reserve down provision percentages from PSI Simulation.
"""
function update_realized_reserve_perc!(device::PSY.Device,
                                        service::PSY.VariableReserve{PSY.ReserveDown},
                                        results::Dict{Symbol, DataFrames.DataFrame},
                                        reserve_up_perc::Dict{String, Array{Float64, 2}},
                                        reserve_down_perc::Dict{String, Array{Float64, 2}})

    reserve_provision = results[Symbol("$(PSY.get_name(service))__VariableReserve_PowerSystems.ReserveDown")][:, Symbol(get_name(device))]
    reserve_perc = reserve_provision / get_device_size(device)
    reserve_down_perc[get_name(device)][1, :] = reserve_perc
    return
end

"""
This function creates the Unit Commitment template for PSI Simulation.
"""
function create_uc_template()
    service = Dict(
        :ReserveUp =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveUp}, PSI.RangeReserve),
        :ReserveDown =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveDown}, PSI.RangeReserve),
    )
    devices = Dict(
        :Generators => PSI.DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment),
        :Ren => PSI.DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch),
        :Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad),
        :HydroROR => PSI.DeviceModel(PSY.HydroEnergyReservoir, PSI.HydroDispatchRunOfRiver),
        :RenFx => PSI.DeviceModel(PSY.RenewableFix, PSI.FixedOutput),
        :Batt => PSI.DeviceModel(PSY.GenericBattery, PSI.BookKeepingwReservation),
    )
    template_uc =
        PSI.template_unit_commitment(network = PSI.PM.NFAPowerModel, devices = devices, services = service)
    return template_uc
end

"""
This function creates the Economic Dispatch template for PSI Simulation.
"""
function create_ed_template()
    service = Dict(
        :ReserveUp =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveUp}, PSI.RangeReserve),
        :ReserveDown =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveDown}, PSI.RangeReserve),
    )
    devices = Dict(
        :Generators => PSI.DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatch),
        :Ren => PSI.DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch),
        :Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad),
        :HydroROR => PSI.DeviceModel(PSY.HydroEnergyReservoir, PSI.HydroDispatchRunOfRiver),
        :RenFx => PSI.DeviceModel(PSY.RenewableFix, PSI.FixedOutput),
        :Batt => PSI.DeviceModel(PSY.GenericBattery, PSI.BookKeepingwReservation),
    )
    template_ed =
        PSI.template_economic_dispatch(network = PSI.PM.NFAPowerModel, devices = devices, services = service)

    return template_ed
end

"""
This function creates the Stages for PSI Simulation.
"""
function create_stage(template::PSI.OperationsProblemTemplate, sys::PSY.System, solver::JuMP.MOI.OptimizerWithAttributes)
    duals = [:nodal_balance_active__Bus]

    services = get_system_services(sys)
    if in(PSY.VariableReserve{PSY.ReserveUp}, typeof.(services))
        push!(duals, Symbol("requirement__VariableReserve_PowerSystems.ReserveUp"))
    end
    if in(PSY.VariableReserve{PSY.ReserveDown}, typeof.(services))
        push!(duals, Symbol("requirement__VariableReserve_PowerSystems.ReserveDown"))
    end

    stage = PSI.Stage(
        PSI.GenericOpProblem,
        template,
        sys,
        solver;
        constraint_duals = duals)

    return stage
end

"""
This function creates the Sequences for PSI Simulation.
"""
function create_sequence(stage_1::AbstractString, stage_2::AbstractString)
    sequence = PSI.SimulationSequence(
        step_resolution = Dates.Hour(24),
        order = Dict(1 => "UC", 2 => "ED"),
        feedforward_chronologies = Dict(("UC" => "ED") => PSI.Synchronize(periods = 24)),
        horizons = Dict("UC" => 24, "ED" => 1),
        intervals = Dict(
            "UC" => (Dates.Hour(24), PSI.Consecutive()),
            "ED" => (Dates.Hour(1), PSI.Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => PSI.SemiContinuousFF(
                binary_source_stage = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = PSI.InterStageChronology(),
    )
    return sequence
end

"""
This function creates the PSI Simulation and post-processes the results.
"""
function create_simulation( sys_UC::PSY.System,
                            sys_ED::PSY.System,
                            zones::Vector{String},
                            days::Int64,
                            solver::JuMP.MOI.OptimizerWithAttributes;
                            kwargs...)

    template_uc = create_uc_template()
    template_ed = create_ed_template()

    uc_stage = create_stage(template_uc, sys_UC, solver)
    ed_stage = create_stage(template_ed, sys_UC, solver)

    stages_definition = Dict(
        "UC" => uc_stage,
        "ED" => ed_stage,
    )

    sequence = create_sequence("UC", "ED")

    sim = PSI.Simulation(
        name = "aggregation",
        steps = days,
        stages = stages_definition,
        stages_sequence = sequence,
        simulation_folder = "./",
    )
    PSI.build!(sim)

    sim_results = PSI.execute!(sim; );

    base_power = PSY.get_base_power(sys_UC)

    res_uc = PSI.load_simulation_results(sim_results, "UC")
    res_ed = PSI.load_simulation_results(sim_results, "ED")

    dual_values = res_ed.dual_values

    duals = keys(dual_values)
    data_length = DataFrames.nrow(dual_values[:dual_nodal_balance_active__Bus])

    energy_price = AxisArrays.AxisArray(zeros(length(zones), 1, data_length), zones, 1:1, 1:data_length)
    reserve_up_price = AxisArrays.AxisArray(zeros(length(zones), 1, data_length), zones, 1:1, 1:data_length)
    reserve_down_price = AxisArrays.AxisArray(zeros(length(zones), 1, data_length), zones, 1:1, 1:data_length)

    for zone in zones
            bus = find_zonal_bus(zone, sys_UC)
            zone_num = parse(Int64, last(zone, 1))
            if isnothing(bus)
                energy_price[zone, 1, :] = zeros(data_length)
            else
                energy_price[zone, 1, :] =
                abs.(round.(dual_values[:dual_nodal_balance_active__Bus][:, PSY.get_name(bus)], digits = 5)) / base_power
            end
            if in(Symbol("dual_requirement__VariableReserve_PowerSystems.ReserveUp"), keys(dual_values))
                reserve_up_price[zone, 1, :] =
                abs.(round.(dual_values[Symbol("dual_requirement__VariableReserve_PowerSystems.ReserveUp")][:, Symbol("Reg_Up_R$(zone_num)")], digits = 5)) / base_power
            end
            if in(Symbol("dual_requirement__VariableReserve_PowerSystems.ReserveDown"), keys(dual_values))
                reserve_down_price[zone, 1, :] =
                abs.(round.(dual_values[Symbol("dual_requirement__VariableReserve_PowerSystems.ReserveDown")][:, Symbol("Reg_Down_R$(zone_num)")], digits = 5)) / base_power
            end
    end

    sys_techs = get_all_techs(sys_UC)
    tech_names = get_name.(sys_techs)
    capacity_factors = Dict([g => zeros(1, data_length) for g in tech_names])
    reserve_up_perc = Dict([g => zeros((1, data_length)) for g in tech_names])
    reserve_down_perc = Dict([g => zeros((1, data_length)) for g in tech_names])

    for tech in sys_techs
        name = get_name(tech)
        services = PSY.get_services(tech)
        capacity_factors[name][1, :] = get_realized_capacity_factors(tech, res_ed.variable_values)

        for service in services
            update_realized_reserve_perc!(tech,
                                         service,
                                         res_ed.variable_values,
                                         reserve_up_perc,
                                         reserve_down_perc)
        end
    end

    return energy_price, reserve_up_price, reserve_down_price, capacity_factors, reserve_up_perc, reserve_down_perc;
end


text:
#SignalProcessing base types & core functions
#-------------------------------------------------------------------------------

#==Useful constants
===============================================================================#
const TIME = DS{:time}()
const FREQ = DS{:freq}()

#==Main data structures
===============================================================================#

#-------------------------------------------------------------------------------
#==TODO: Remove.  Keep around in case is still useful.
abstract type AbstractDataDomain end
abstract type Domain{Symbol} <: AbstractDataDomain end

#Aliases to reduce typo mistakes:
const BitDomain = Domain{:bit} #Bit domain (discrete time, but one sample per bit)
#TODO: Should BitDomain <: DTDomain?  Should it not exist?
const DTDomain = Domain{:DT} #Discrete-time domain
const CTDomain = Domain{:CT} #Continuous-time domain
const FDomain = Domain{:f} #Frequency domain
const ZDomain = Domain{:Z} #Z-domain
==#

#-------------------------------------------------------------------------------
abstract type AbastractFrequency{Symbol} end #Id
const Hertz = AbastractFrequency{:Hz}
const RadiansPerSecond = AbastractFrequency{:rad}

#-------------------------------------------------------------------------------
struct Pole{T<:Number,AbastractFrequency}
	v::T
end
const PoleHz{T<:Number} = Pole{T,Hertz}
const PoleRad{T<:Number} = Pole{T,RadiansPerSecond}

#Constructor function:
Pole(v::T, u::Symbol) where {T<:Number} = Pole{T,AbastractFrequency{u}}(v)

#==Time<=>Frequency domain signals
Only supports real signals for now (using rfft).
Users need only provide positive frequency spectrum.
TODO: Could we make this more generic (rename?) to include pos<=>"Frequency"?
TODO: What about 2D images, 3D movies, etc?
TODO: Make generic, instead of basing t/f on DataFloat?
TODO: Is DataTF necessary??  Could having just DataTime & DataFreq be enough?
===============================================================================#
mutable struct DataTF
	tperiodic::Bool #Represents a periodic signal
	tvalid::Bool #Time-domain data is valid
	fvalid::Bool #Frequency-domain data is valid
	#NOTE: Keep length(t) as either odd/even, when performing rfft
	t::StepRangeLen{DataFloat}
	xt::Vector{DataFloat} #x(t)
	f::StepRangeLen{DataFloat}
	Xf::Vector{DataComplex} #X(f)
end
DataTF(::DS{:time}, t, xt; tperiodic::Bool=false) =
	DataTF(tperiodic, true, false, t, xt, timetofreq(t), [])

DataTF(::DS{:freq}, f, Xf; teven::Bool=true, tperiodic::Bool=false) =
	DataTF(tperiodic, false, true, freqtotime(f, teven), [], f, Xf)

#Simple wrapper to point to DataTF:
#-------------------------------------------------------------------------------
mutable struct DataTime
	data::DataTF
end
DataTime(t::AbstractRange, xt::Vector{DataFloat}; tperiodic::Bool=false) = DataTime(DataTF(TIME, t, xt, tperiodic=tperiodic))
DataTime(t::AbstractRange; tperiodic::Bool=false) = DataTime(t, collect(t), tperiodic=tperiodic)
#For low-level algorithms... Initializes data ranges:
DataTime(d::DataTF, xt::Vector{DataFloat}) =
	validatelengths(DataTime(DataTF(d.tperiodic, true, false, d.t, xt, d.f, [])))
DataTime(t::AbstractRange, xt::DataF1; tperiodic::Bool=false) =
	DataTime(t, sample(xt, t).y, tperiodic=tperiodic)

mutable struct DataFreq
	data::DataTF
end
DataFreq(f::AbstractRange, Xf::Vector{DataComplex}; teven::Bool=true, tperiodic::Bool=false) =
	DataFreq(DataTF(FREQ, f, Xf, teven=teven, tperiodic=tperiodic))
DataFreq(f::AbstractRange; teven::Bool=true, tperiodic::Bool=false) =
	DataFreq(f, collect(f), teven=teven, tperiodic=tperiodic)
DataFreq(f::AbstractRange, Xf::DataF1; teven::Bool=true, tperiodic::Bool=false) =
	DataFreq(f, sample(Xf, f).y, teven=teven, tperiodic=tperiodic)
#For low-level algorithms... Initializes data ranges:
DataFreq(d::DataTF, Xf::Vector{DataComplex}) =
	validatelengths(DataFreq(DataTF(d.tperiodic, false, true, d.t, [], d.f, Xf)))

#Z-domain data (discrete time, continuous frequency approximation):
#TODO: Think about details a bit more... what about x-values, etc?
DataZ(n::Int; teven::Bool=true) =
	DataFreq(linspace(0, 2pi, n), Xf, teven=teven, tperiodic=false)
DataZ(Xf::Vector{DataComplex}; teven::Bool=true) =
	DataFreq(linspace(0, 2pi, length(Xf)), Xf, teven=teven, tperiodic=false)

#==Useful validations/assertions
===============================================================================#

#Validate data lengths:
function validatelengths(d::DataTime)
	ensure(length(d.data.t)==length(d.data.xt),
		ArgumentError("Invalid DataTime: Lengths do not match: t & x(t)"))
	return d
end
function validatelengths(d::DataFreq)
	ensure(length(d.data.f)==length(d.data.Xf),
		ArgumentError("Invalid DataFreq: Lengths do not match: f & X(f)"))
	return d
end

#==Conversion functions
===============================================================================#
#==Problem:
Conversions do not maintain numeric data types.
TODO: Use proper type conversion routines
==#
Base.convert(::Type{U}, p::Pole{T,U}) where {T,U<:AbastractFrequency} = p
Base.convert(::Type{Hertz}, p::PoleRad) = Pole(p.v/(2*pi),:Hz)
Base.convert(::Type{RadiansPerSecond}, p::PoleHz) = Pole(p.v*(2*pi),:rad)

#Value extractor functions:
value(p::Pole) = p.v
value(u::Symbol, p::Pole) = value(convert(AbastractFrequency{u}, p))

#==Main functions
===============================================================================#

#Last line


text:
const TIMELINE_SEQUENCE = 1
const TIMELINE_PARALLEL = 2

# Used as parameter in [repeat] and [repeatYoyo] methods.
const INFINITY = -1

# A typical max would be 4, for example a color would need an array size
# of [r, g, b, a] = 4.
# a position [x,y] or alpha [a] would certainly fit.
const DEFAULT_COMBINED_ATTRS_LIMIT = 4
const EPSILON = 0.00000000001

tween_pool = TweenPool()
accessors = Dict{Int64,AbstractTweenAccessor}()

mutable struct Tween <: AbstractTween
    base::BaseTween

    easing::AbstractTweenEquation

    # An object associated with this tween.
    # The target is what is read/write via the accessor
    target::Any

    # Class types could be: 1 = AbstractNode, 2 = Vector2D, 3 = Color etc...
    class_type::Int64 # (aka tween_type)

    # Each accessor handles a different property on the target.
    accessor::AbstractTweenAccessor

    # Indicates how many values are to be modified by the currrent update time step.
    combined_attrs_cnt::Int64
    combined_attrs_limit::Int64

    # Indicates which value we want the accessor to read or write.
    value_type::Int64

    start_values::Array{Float64,1}
    # These are the values the tween interpolates to (aka targets for).
    # It is set by the user while creating the tween. For example, to tween from
    # black (0,0,0) to red (255,0,0) the target values would be:
    # [255, 0, 0] array.
    # Code example:
    #   tween.target_values[1] = 255
    #   tween.target_values[2] = 0
    #   tween.target_values[3] = 0
    target_values::Array{Float64,1}
    accessor_buffer::Array{Float64,1}

    # General
    is_from::Bool
    is_relative::Bool

    function Tween()
        o = new()
        o.base = BaseTween()

        o.combined_attrs_limit = DEFAULT_COMBINED_ATTRS_LIMIT

        o.start_values = Array{Float64}(undef, DEFAULT_COMBINED_ATTRS_LIMIT)
        fill!(o.start_values, 0.0)

        o.target_values = Array{Float64}(undef, DEFAULT_COMBINED_ATTRS_LIMIT)
        fill!(o.target_values, 0.0)

        o.accessor_buffer = Array{Float64}(undef, DEFAULT_COMBINED_ATTRS_LIMIT)
        fill!(o.accessor_buffer, 0.0)

        o.target = nothing
        o
    end
end

function pool_size()
    count(tween_pool)
end

# Factories -------------------------------------------------------------
# value_type = tween_type
function create_tween_to(node::Ranger.AbstractNode, value_type::Int64, duration::Float64)
    tween = get_from!(tween_pool)

    # Be sure to change this default easing to match your requirements.
    tween.easing = Linear(EASE_INOUT) # A default equation

    setup!(tween, node, value_type, duration)

    tween
end
# ------------------------------------------------------------------------

function register_accessor(class_type::Int64, accessor::AbstractTweenAccessor)
    accessors[class_type] = accessor
end

# return a AbstractTweenAccessor
function get_registered_accessor(class_type::Int64)
    accessors[class_type]
end

# Changes the [limit] for combined attributes. Defaults to 3 to reduce memory footprint.
function set_combined_attrs_limit(tween::AbstractTween, limit::Int64)
    tween.combined_attrs_limit = limit
end

function start!(tween::AbstractTween)
    if tween.base.auto_start
        build!(tween)
        start!(tween.base)
    end
end

function reset!(tween::AbstractTween)
    reset!(tween.base)
    println("tween resetting")
    tween.target = nothing
    tween.class_type = -1
    tween.value_type = -1
    tween.easing = nothing
    tween.is_from = false
    tween.is_relative = false
    tween.combined_attrs_cnt = 0

    if length(tween.accessor_buffer) ≠ tween.combined_attrs_limit
        tween.accessor_buffer = Array{Float64}(undef, tween.combined_attrs_limit)
        fill!(tween.accessor_buffer, 0.0)
    end
end

function update!(tween::AbstractTween, delta::Float64)
    update!(tween.base, delta)
end

function kill!(tween::AbstractTween)
    kill!(tween.base)
end

function reset!(tween::AbstractTween)
    reset!(tween.base)
end

# Sets the [TweenCallbackHandler]. By default, it will be fired at the completion of the tween or timeline (event COMPLETE).
# If you want to change this behavior and add more triggers, use the [setCallbackTriggers] method.
function set_callback!(tween::AbstractTween, callback::Function)
    tween.base.state_callback = callback
end

function set_target_values!(tween::AbstractTween, values::Array{Float64,1})
    idx = 1
    for value in values
        tween.target_values[idx] = value
        idx += 1
    end
end

function callCallback(tween::AbstractTween, type::UInt64)
    if !isnothing(tween.base.state_callback) && (tween.base.callback_triggers & type) > 0
        tween.base.state_callback(type, tween)
    end
end

function build!(tween::AbstractTween)
    if isnothing(tween.target)
        return
    end

    tween.accessor = get_registered_accessor(tween.class_type)

    if isnothing(tween.accessor)
        tween.accessor = tween.target
    end

    if isnothing(tween.accessor)
        throw("No TweenAccessor was found for the target, and it is not Tweenable either.")
    end

    # Fetch the current range values from the target.
    tween.combined_attrs_cnt =
        get_values!(tween.accessor, tween, tween.accessor_buffer) # get_values will use 'tween.value_type' as a property selector

    if tween.combined_attrs_cnt < 0
        tween.combined_attrs_cnt = 0
    end

    if tween.combined_attrs_cnt > tween.combined_attrs_limit
        msg = """
        You cannot combine more than 'combined_attrs_limit'
        attributes in a tween. You can raise this limit with
        Tween.setCombinedAttributesLimit(), which should be called once
        in application initialization code.
        """
        throw(msg)
    end
end

function initialize!(tween::AbstractTween)
    if tween.base.current_time + tween.base.delta_time >= tween.base.delay
        initialize_override!(tween)
        tween.base.initialized = true
        tween.base.is_iterationStep = true
        tween.base.step = 0
        tween.base.delta_time -= tween.base.delay - tween.base.current_time
        tween.base.current_time = 0
        tween.base.state_callback(TWEEN_CALLBACK_BEGIN, tween)
        tween.base.state_callback(TWEEN_CALLBACK_START, tween)
    end
end

function initialize_override!(tween::AbstractTween)
    if (tween.target == nothing)
        return
    end

    get_tweened_values!(tween, tween.start_values)

    for i in 1:tween.combined_attrs_cnt
        inc = if tween.is_relative
            tween.start_values[i]
        else
            0
        end
        tween.target_values[i] += inc

        if tween.is_from
            tmp = tween.start_values[i]
            tween.start_values[i] = tween.target_values[i]
            tween.target_values[i] = tmp
        end
    end

end

function update_override!(tween::AbstractTween, step::Int64, last_step::Int64, is_iteration_step::Bool, delta::Float64)
    if (tween.target == nothing || tween.easing == nothing)
        return
    end

    # in case iteration end has been reached
    if !is_iteration_step && step > last_step
        values = if is_reverse(tween, last_step)
            tween.start_values
        else
            tween.target_values
        end

        set_tweened_values!(values)
        return
    end

    if !is_iteration_step && step < last_step
        values = if is_reverse(tween, last_step)
            tween.target_values
        else
            tween.start_values
        end

        set_tweened_values!(values)
        return
    end

    # Validation
    @assert is_iteration_step "Expected an iteration step"
    @assert tween.base.current_time >= 0.0 "Expected current time > 0.0"
    @assert tween.base.current_time <= tween.base.duration "Expected an current time < duration"

    # Case duration equals zero
    if tween.base.duration < EPSILON && delta > -EPSILON
        values = if is_reverse(tween, step)
            tween.target_values
        else
            tween.start_values
        end

        set_tweened_values!(values)
        return
    end

    if tween.base.duration < EPSILON && delta < EPSILON
        values = if is_reverse(tween, step)
            tween.start_values
        else
            tween.target_values
        end

        set_tweened_values!(values)
        return
    end

    # Normal behavior
    time = if is_reverse(tween, step)
        tween.base.duration - tween.base.current_time
    else
        tween.base.current_time
    end

    # compute modifies "t" based on easing equations. This value
    # is then passed to the lerp loop below.
    t = tween.easing.compute(time / tween.base.duration)

    for i in 1:tween.combined_attrs_cnt
        tween.accessor_buffer[i] =
            tween.start_values[i] + t * (tween.target_values[i] - tween.start_values[i])
    end

    set_tweened_values!(tween, tween.accessor_buffer)
end

function force_start_values!(tween::AbstractTween)
    if tween.target == nothing
        return
    end
    set_tweened_values!(tween, tween.start_values)
end

function force_end_values!(tween::AbstractTween)
    if tween.target == nothing
        return
    end
    set_tweened_values!(tween, tween.target_values)
end

# class_type = tween_type
function contains_target(tween::AbstractTween, target::AbstractTween, class_type::Int64)
    tween.target == target && tween.class_type == class_type
end

# Stops and resets the tween or timeline, and sends it to its pool, for
# later reuse. Note that if you use a [TweenManager], this method
# is automatically called once the animation is finished.
function free!(tween::AbstractTween)
    put_to!(tween_pool, tween)
end

function setup!(tween::AbstractTween, node::Ranger.AbstractNode, value_type::Int64, duration::Float64)
    if duration < 0.0
        throw("Duration can't be negative")
    end

    tween.value_type = value_type
    tween.base.target = node
    tween.base.duration = duration
end

function get_tweened_values!(tween::AbstractTween, into_buffer::Array{Float64,1})
    value_count = if tween.accessor ≠ nothing
        get_values!(tween.accessor, tween, into_buffer)
    else
        0
    end

    value_count
end

function set_tweened_values!(tween::AbstractTween, values::Array{Float64,1})
    if tween.accessor ≠ nothing
        set_values!(tween.accessor, tween, values)
    end
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
function rodcut(totlength::Int, prices::Array{Tuple{Int, Float}, 1})
	price_list::Array{Float, 1} = zeros(Float, totlength)

	for i in 1 : length(prices)
		price_list[prices[i][1]] = prices[i][2]
	end
	rodcut(totlength, price_list)
	return price_list[totlength]
end

function rodcut(totlength::Int, prices::Array{Float, 1})
	for i in 1 : totlength
		for j in 1 : i - 1
			if prices[j] + prices[i - j] > prices[i]
				prices[i] = prices[j] + prices[i - j]
			end
		end
	end
end

function rodcut(totlength::Float, prices::Array{Tuple{Float, Float}})

end


text:
abstract type AbstractTrainableDist end

"""
    TrainableMvNormal(shape;
                    init=glorot_normal,
                    device=cpu) <: AbstractTrainableDist
    TrainableMvNormal(mean, stddev, sample, shape)

A Multivariate Normal distribution with trainable mean and stddev.

# Fields
- `mean`: Trainable mean vector of the distribution
- `stddev`: Trainable standard deviation vector of the distibution
- `sample`: The latest sample from the distribution, used in calculating loglikelhood loss
- `shape::Tuple`: The shape of the sample to be returned

# Arguments
- `shape::Tuple`: The shape of the sample to returned from the distribution
- `init`: glorot_normal; to initialize the mean and stddev trainable params
- `device`: cpu; the device to move the sample to, used for convinience while using both GPU and CPU
"""
struct TrainableMvNormal{M,S,N,D} <: AbstractTrainableDist
    mean::M
    stddev::M
    sample::S
    shape::NTuple{N,Integer}
    device::D
end

function TrainableMvNormal(shape; init = glorot_normal, device = cpu)
    # Create the mean and stddev
    total_params = prod(shape)
    mean = device(init(total_params))
    stddev = device(init(total_params))
    sample = device(zeros(total_params))
    return TrainableMvNormal(mean, stddev, sample, shape, device)
end

# Don't backprop through the sample stored for loss calc
@functor TrainableMvNormal (mean, stddev)

"""
    Forward pass of the trainable distribution, returns a sample from
a multivariate normal with the trained mean and stddev
"""
function (tmv::TrainableMvNormal)()
    dist = DistributionsAD.TuringMvNormal(tmv.mean, abs.(tmv.stddev))
    # Sample from the dist
    # Put it on zygote.ignore to supress mutation errors
    sample = tmv.device(rand(dist))
    # Ignore the mutation array error while backprop
    Flux.Zygote.@ignore copyto!(tmv.sample, sample)
    return reshape(sample, tmv.shape)
end


text:
#===============================================================================
# Practice implementing a classic data structure in Julia
#
# Author: Carl Cortright
# Date: 1/1/2017
#
===============================================================================#

__precompile__()

module linkedlist

export DoubleLinkedList, append, insert, remove, listsize, access

#
# Each node in the linked list represents some data
#
type Node{T}
  next::Node{T}
  prev::Node{T}
  data::T

  Node() = (self = new(); self.next = self; self.prev = self)
end

#
# The main linked list object
#
type DoubleLinkedList{T}
  head::Node{T}
  tail::Node{T}
  size::Int

  function DoubleLinkedList()
    head = Node{T}()
    tail = head
    head.next = tail
    new(head, tail, 0)
  end
end

#
# Adds a node with data to the end of the list
#
function append{T}(list::DoubleLinkedList{T}, data::T)::Node{T}
  newNode = Node{T}()
  newNode.data = data
  newNode.prev = list.tail
  newNode.next = list.head
  list.head.prev = newNode
  list.tail.next = newNode
  list.tail = newNode
  list.size += 1
  return newNode
end

#
# Adds a node with data at position
#
function insert{T}(list::DoubleLinkedList{T}, data::T, position::Int)::Node{T}
  if position > list.size
    error("Attempting to access position $position which is out of range.")
  end
  current = list.head
  for i in 1:position
    current = current.next
  end
  insertNode = Node{T}()
  insertNode.data = data
  insertNode.prev = current.prev
  insertNode.next = current
  current.prev.next = insertNode
  current.prev = insertNode
  list.size += 1
  return insertNode
end

#
# Removes the node at position
#
function remove{T}(list::DoubleLinkedList{T}, position::Int)::T
  if position > list.size
    error("Attempting to access position $position which is out of range.")
  end
  current = list.head
  for i in 1:position
    current = current.next
  end
  data = current.data
  current.prev.next = current.next
  current.next.prev = current.prev
  list.size -= 1
  return data
end

#
# Gets the data from the node at position
#
function access{T}(list::DoubleLinkedList{T}, position::Int)::T
  if position > list.size
    error("Attempting to access position $position which is out of range.")
  end
  current = list.head
  for i in 1:position
    current = current.next
  end
  return current.data
end

#
# Gets the size of the list
#
function listsize{T}(list::DoubleLinkedList{T})
  return list.size
end

end  # module LinkedList

max_issues_repo_name: carlcortright/Playground


text:
module Dummy1

using ComputationalResources

function __init__()
    push!(LOAD_PATH, dirname(@__FILE__))
    if haveresource(ArrayFireLibs)
        @eval using Dummy1AF
    end
    pop!(LOAD_PATH)
end

foo(::CPU1) = true

end


text:
using ClangCompiler
using ClangCompiler.LLVM
using Test

@testset "JIT call" begin
    src = joinpath(@__DIR__, "code", "main.cpp")
    args = get_compiler_args()
    haskey(ENV, "CI") && push!(args, "-v")
    irgen = IRGenerator(src, args)
    ee = JIT(get_module(irgen))
    f = lookup_function(ee, "main")
    ret = convert(Cint, run(ee, f))
    @test ret == 42
    dispose(irgen)
end


text:
##

function _locpoly_alt1!(
  eX::Matrix{T}, w::Vector{T},
  X::Matrix{T}, z::Vector{T}, y::Vector{T},
  z0::T, degree::Int64, kernel::SmoothingKernel{T}) where {T <: AbstractFloat}

  @. w = evaluate(kernel, z, z0)
  _expand_X!(eX, X, z, z0, degree)

  ((eX' * diagm(w)) * eX) \ ((eX' * diagm(w)) * y)
end

function locpoly_alt1(
  X::Matrix{T}, z::Vector{T}, y::Vector{T},
  zgrid::Vector{T},
  degree::Int64,                                     # degree of the polynomial
  kernel::SmoothingKernel{T}=GaussianKernel(one(T))) where {T <: AbstractFloat}

  n, p = size(X)
  ep = p * (degree + 1)
  out = Array{T}(ep, length(zgrid))
  w = Array{T}(n)
  eX = Array{T}(n, ep)

  ind = 0
  for z0 in zgrid
    ind += 1
    out[:, ind] = _locpoly_alt1!(eX, w, X, z, y, z0, degree, kernel)
  end
  out
end

##

function _locpoly_alt!(
  Xt_w_X::Matrix{T}, Xt_w_Y::Vector{T}, w::Vector{T},
  X::Matrix{T}, z::Vector{T}, y::Vector{T},
  z0::T, degree::Int64, kernel::SmoothingKernel{T}) where {T <: AbstractFloat}

  @. w = evaluate(kernel, z, z0)
  _expand_Xt_w_X!(Xt_w_X, w, X, z, z0, degree)
  _expand_Xt_w_Y!(Xt_w_Y, w, X, z, y, z0, degree)

  cholfact!(Xt_w_X) \ Xt_w_Y
end

function locpoly_alt(
  X::Matrix{T}, z::Vector{T}, y::Vector{T},
  zgrid::Vector{T},
  degree::Int64,                                     # degree of the polynomial
  kernel::SmoothingKernel{T}) where {T <: AbstractFloat}

  n, p = size(X)
  ep = p * (degree + 1)
  out = Array{T}(ep, length(zgrid))
  w = Array{T}(n)
  Xt_w_X = Array{T}(ep, ep)
  Xt_w_Y = Array{T}(ep)

  ind = 0
  for z0 in zgrid
    ind += 1
    out[:, ind] = _locpoly_alt!(Xt_w_X, Xt_w_Y, w, X, z, y, z0, degree, kernel)
  end
  out
end

##

function locpolyl1_alt(
  X::Matrix{T}, z::Vector{T}, y::Vector{T},
  zgrid::Vector{T},
  degree::Int64,
  kernel::SmoothingKernel{T},
  λ0::T,
  options::CDOptions=CDOptions()) where {T <: AbstractFloat}

  # construct inner options because we do not want warmStart = false
  # we want to start from the previous iteration since the points
  # on the grid should be close to each other
  opt = CDOptions(options.maxIter, options.optTol, options.randomize, true, options.numSteps)

  n, p = size(X)
  ep = p * (degree + 1)
  out = Array{SparseVector{T, Int64}}(length(zgrid))

  # temporary storage
  w = Array{T}(n)
  wX = Array{T}(n, ep)
  stdX = ones(T, ep)

  ind = 0
  for z0 in zgrid
    ind += 1

    # the following two should update f
    @. w = evaluate(kernel, z, z0)
    _expand_X!(wX, X, z, z0, degree)
    # compute std for each column
    @inbounds for j=1:ep
      v = zero(T)
      @simd for i=1:n
        v += wX[i, j] * wX[i, j] * w[i]
      end
      stdX[j] = sqrt(v / n)
    end

    @. w = sqrt(w)
    f = CDLeastSquaresLoss(diagm(w)*y, diagm(w)*wX)
    g = ProximalBase.ProxL1(λ0, stdX)
    β = ProximalBase.SparseIterate(ep)

    # solve for β
    coordinateDescent!(β, f, g, opt)
    out[ind] = convert(SparseVector, β)
  end
  out
end

##

function _expand_Xt_w_X_alt!(
  Xt_w_X::Matrix{T},
  w::Vector{T}, X::Matrix{T},
  z::Vector{T}, z0::T, degree::Int64) where {T <: AbstractFloat}

  n, p = size(X)
  fill!(Xt_w_X, zero(T))
  @inbounds for j=1:p, jj=0:degree
    col=(j-1)*(degree+1)+1+jj
    for k=j:p
      if k != j
        krange = 0:degree
      else
        krange = jj:degree
      end
      for kk=krange
        row = (k-1)*(degree+1)+1+kk

        # compute Xt_w_X[row, col]
        for i=1:n
          Xt_w_X[row, col] += X[i, j] * (z[i] - z0)^(jj+kk) * X[i, k] * w[i]
        end
        if row != col
          Xt_w_X[col, row] = Xt_w_X[row, col]
        end
      end
    end
  end
  Xt_w_X
end

function genData(n, p)
  X = randn(n, p)
  Z = rand(n)
  ɛ = 0.1 * randn(n)
  Y = zeros(n)
  betaMat = zeros(p, n)

  randBeta = [rand([2,4,6,8]) for j=1:p]
  for i=1:n
    betaMat[:, i] = sin.(randBeta * Z[i])
    Y[i] = dot(betaMat[:, i], X[i, :]) + ɛ[i]
  end
  Y, X, Z, betaMat
end

function genDataIID(n, p)
  X = randn(n, p)
  Z = rand(n)
  ɛ = 0.1 * randn(n)
  Y = zeros(n)
  beta = 0.5+rand(p)

  for i=1:n
    Y[i] = dot(beta, X[i, :]) + ɛ[i]
  end
  Y, X, Z, beta
end

# using BenchmarkTools

####
#
# n, p = 100, 2
# Y, X, Z, betaMat = genData(n, p)
# zgrid = collect(0.01:0.1:0.99)
#
# gk = GaussianKernel(0.2)
# degree = 3
#
# @time o1 = locpoly(X, Z, Y, zgrid, 1, gk)
# @time o2 = locpoly_alt(X, Z, Y, zgrid, 1, gk)
# @time o3 = locpoly_alt1(X, Z, Y, zgrid, 1, gk)
#
# maximum(abs.(o1-o2))
# maximum(abs.(o1-o3))
#
# @benchmark locpoly($X, $Z, $Y, $zgrid, $degree, $gk)
# @benchmark locpoly_alt($X, $Z, $Y, $zgrid, $degree, $gk)
# @benchmark locpoly_alt1($X, $Z, $Y, $zgrid, $degree, $gk)

####

# p = 10
# X = randn(100, p)
# z = rand(100)
# w = zeros(100)
# k = GaussianKernel(0.2)
# @. w = evaluate(k, z, 0.5)
#
#
# degree = 2
# cp = p*(degree+1)
#
# eX = zeros(100, cp)
# _expand_X!(eX, X, z, 0.5, degree)
#
# Xt_w_X = zeros(cp, cp)
# Xt_w_X1 = zeros(cp, cp)
# o1 = _expand_Xt_w_X!(Xt_w_X, w, X, z, 0.5, degree)
# o2 = _expand_Xt_w_X_alt!(Xt_w_X1, w, X, z, 0.5, degree)
# o3 = (eX'*diagm(w))*eX
#
# maximum(abs.(o1 - o2))
# maximum(abs.(o1 - o3))
#
# using BenchmarkTools
# @benchmark _expand_Xt_w_X!($Xt_w_X, $w, $X, $z, 0.5, degree)
# @benchmark _expand_Xt_w_X_alt!($Xt_w_X1, $w, $X, $z, 0.5, degree)

### benchmark local lasso
#
#
# n, s, p = 500, 2, 50
# degree = 2
# gk = GaussianKernel(0.2)
# w = zeros(n)
# sw = zeros(n)
#
#
# Y, X, Z, betaMat = genData(n, s)
# X = [X zeros(n, p-s)]
# z0 = rand()
#
# cp = p*(degree+1)
# eX = zeros(n, cp)
# CoordinateDescent._expand_X!(eX, X, Z, z0, degree)
#
# @. w = evaluate(gk, Z, z0)
# @. sw = sqrt(w)
#
# λ = 0.1
# g = ProximalBase.ProxL1(λ)
# f1 = CDLeastSquaresLoss(diagm(sw) * Y, diagm(sw)*X)
# f2 = CDWeightedLSLoss(Y, X, w)
#
# opt = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=false, randomize=true)
# x1 = convert(ProximalBase.SparseIterate, sprand(p, 0.6))
# x2 = convert(ProximalBase.SparseIterate, sprand(p, 0.6))
#
# @benchmark coordinateDescent!($x1, $f1, $g, $opt)
# @benchmark coordinateDescent!($x2, $f2, $g, $opt)

## benchmark locpolyl1
#
# n, s, p = 500, 2, 50
# degree = 0
# gk = GaussianKernel(0.2)
#
# Y, X, Z, betaT = genData(n, s)
# X = [X zeros(n, p-s)]
# zgrid = collect(0.01:0.1:0.99)
#
# λ0 = 0.01
#
# opt = CDOptions(;randomize=false)
# o1 = locpolyl1(X,Z,Y,zgrid,degree,gk,λ0, opt)
# o2 = locpolyl1_alt(X,Z,Y,zgrid,degree,gk,λ0, opt)
#
# maximum( maximum(abs.(o1[i] - o2[i])) for i=1:length(zgrid) )

#


text:
function initialize_outer!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicInverter{
        C,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
        IC,
        DC,
        P,
        F,
    },
) where {
    C <: PSY.Converter,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(
        dynamic_device,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
    )
    Vr_filter = device_states[external_ix[4]]
    Vi_filter = device_states[external_ix[5]]
    Ir_filter = device_states[external_ix[6]]
    Ii_filter = device_states[external_ix[7]]

    Vr_cnv = get_inner_vars(dynamic_device)[Vr_cnv_var]
    Vi_cnv = get_inner_vars(dynamic_device)[Vi_cnv_var]
    θ0_oc = angle(Vr_cnv + 1im * Vi_cnv)

    #Obtain additional expressions
    p_elec_out = Ir_filter * Vr_filter + Ii_filter * Vi_filter
    q_elec_out = -Ii_filter * Vr_filter + Ir_filter * Vi_filter

    #Update inner_vars
    get_inner_vars(dynamic_device)[P_ES_var] = p_elec_out
    #Update states
    outer_ix = get_local_state_ix(
        dynamic_device,
        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
    )
    outer_states = @view device_states[outer_ix]
    outer_states[1] = PSY.get_ω_ref(dynamic_device) #ω
    outer_states[2] = θ0_oc #θ_oc
    outer_states[3] = q_elec_out #qm

    #Update inner vars
    get_inner_vars(dynamic_device)[θ_oc_var] = θ0_oc
    get_inner_vars(dynamic_device)[ω_oc_var] = PSY.get_ω_ref(dynamic_device)
    #Update Q_ref. Initialization assumes q_ref = q_elec_out of PF solution
    PSY.get_ext(dynamic_device)[CONTROL_REFS][P_ref_index] = p_elec_out
    PSY.set_P_ref!(PSY.get_active_power(PSY.get_outer_control(dynamic_device)), p_elec_out)
    PSY.get_ext(dynamic_device)[CONTROL_REFS][Q_ref_index] = q_elec_out
end


text:
function advancetime!(system::LegSystem,spl::Dierckx.Spline1D,times::CVTimer,
    n::Int64,nsamp::Int64,terms::Vector{Int64},splits::Vector{Int64};runtype="predict")
    tic();
    if runtype == "predict"
        while system.solverparams.numbeats < system.solverparams.numbeatstotal
            # time within heart cycle
            tp = system.t[n+1] - sum(system.solverparams.th*system.solverparams.numbeats);
            # interpolate proximal flow velocity from patient data
            uprox = spl(tp) #+ rand(Distributions.Normal(0,0.05));
            # TVD RK3 time integration
            LegHemodynamics.tvdrk3!(system,times,n,splits,terms,-uprox);
            tic();
            # arterial pressure update
            LegHemodynamics.arterialpressure!(system,n);
            # check for start of new cardiac cycle
            oldnumbeats = system.solverparams.numbeats;
            LegHemodynamics.setnumbeats!(system,n);
            if oldnumbeats < system.solverparams.numbeats
                # update time discretization
                if system.solverparams.numbeats < system.solverparams.numbeatstotal
                    LegHemodynamics.updatediscretization!(system);
                end
            end
            times.tr += toq();
            if mod(n/nsamp,1) == 0
                n+=1;
                break
            end
            n+=1;
        end
    elseif runtype == "assim"
        # check if new cardiac cycle has just begun
        numbeats = system.solverparams.numbeats;
        ttotal = system.solverparams.th*(system.solverparams.numbeats+1);
        if system.t[n+1] > ttotal
            println("New cardiac cycle beginning. Decrementing number of cycles for assimilation.")
            numbeats -= 1;
        end
        # time within heart cycle
        tp = system.t[n] - sum(system.solverparams.th*numbeats);
        # interpolate proximal flow velocity from patient data
        uprox = spl(tp) #+ rand(Distributions.Normal(0,0.05));
        # TVD RK3 time integration
        LegHemodynamics.tvdrk3!(system,times,n-1,splits,terms,-uprox);
        tic();
        # arterial pressure update
        LegHemodynamics.arterialpressure!(system,n-1);
        times.tr += toq();
    end
    if runtype == "predict"
        times.tt += toq();
    elseif runtype == "assim"
        times.tt = toq();
    end
    return system,n,times
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
struct Protected end

protect(x) = Sym{Protected}(:protected, x)

unprotect(x) = x
unprotect(x::Sym{TAG}) where {TAG} = Sym{TAG}(gethead(x), map(unprotect, getargs(x))...)
unprotect(x::Sym{Protected}) = unprotect(firstarg(x))

isprotected(x) = false
isprotected(x::Sym{Protected}) = true


text:
module TestRemat

using Base.Test, MixedModels, FMRIModels, DataFrames

# from MixedModels/test/data.jl:
const ds = DataFrame(Yield = [1545.,1440.,1440.,1520.,1580.,1540.,1555.,1490.,1560.,1495.,
                              1595.,1550.,1605.,1510.,1560.,1445.,1440.,1595.,1465.,1545.,
                              1595.,1630.,1515.,1635.,1625.,1520.,1455.,1450.,1480.,1445.],
                     Batch = pool(repeat('A' : 'F', inner = 5)))

m = lmm(@formula(Yield ~ 1 + (1|Batch)), ds)

terms = m.trms

erm = ExpandedReMat(terms[1])

@test sparse(erm) == sparse(terms[1])

fit!(m)

m2 = lmm(@formula(Yield ~ 1 + (1|Batch)), ds)
m2.trms[1] = erm

fit!(m2)

end # module


text:
module TestExamplesTutorialMissings
using LiterateTest.AssertAsTest: @assert
include("../examples/tutorial_missings.jl")
end  # module


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbeval(nurbs::NURBS{I,F},
            tt::Union{Vector{F}, Vector{Vector{F}}, Matrix{F}} [,
            flag::Symbol=:cartesian]
           ) where {I<:Integer, F<:AbstractFloat}

Evaluate a NURBS structure at parametric points.

Evaluation of NURBS curves or surfaces at parametric points along the U and V
directions. Either homogeneous coordinates are returned if the flag is set
accordingly otherwise cartesian coordinates. This function utilises the
function [`bspeval`](@ref).

# Arguments:
- `nurbs`: NURBS structure
- `tt`: parametric evaluation points either as `Vector{F}` or
  `Vector{Vector{F}}`
- `flag`: flag indicating that cartesian coordinates or homogeneous coordinates
  should be returned (valid values: `:cartesian` (default) or `:homogeneous`)

# Output:
- `p::Union{Matrix{F},Array{F,3}}`: Evaluated points on the NURBS curve or
   surface as cartesian coordinates (x,y,z). If `flag=:homogeneous` is set the
   points are returned as homogeneous coordinates (wx,wy,wz).
- `w::Union{Vector{F},Matrix{F}}`: Weights of the homogeneous coordinates of
  the evaluated points. Note: Only returned it the corresponding flag is set
  (see above).

# Examples:
```julia
julia> p = nrbeval(nurbs,tt)
julia> p,w = nrbeval(nurbs,tt,:homogeneous)
julia> p = nrbeval(srf,[ut,vt])
```

Evaluate the NURBS circle with twenty points from 0.0 to 1.0.
```julia
julia> nrb = nrbcirc();
julia> ut = collect(range(0.0,stop=1.0,length=20))
julia> p = nrbeval(nrb,ut)
```

!!! note "Note:"
    Contrary to the Matlab version does this ported version evaluates not the
    `varargout` count to switch between different return values. A symbol
    `flag` is used instead.

# Reference
- see also: [`bspeval`](@ref)
"""
function nrbeval(nurbs::NURBS{I,F},
                 tt::Union{Vector{F}, Vector{Vector{F}}, Matrix{F}},
                 flag::Symbol=:cartesian
                ) where {I<:Integer,F<:AbstractFloat}
# TODO: this ported version does not do all the tests of the matlab version,
#       instead the type features of julia are used -> check if the matlab
#       errors get caught

# check if flag is a valid symbol
if !(flag in [:cartesian;:homogeneous])
  throw(ArgumentError("Valid flags are :cartesian or :homogeneous !"));
end # if

if typeof(nurbs)<:NURBS2D
  if typeof(tt)<:Vector{F}
    throw(ArgumentError("Parametric evaluation points for surfaces are " *
                        "either Vector{Vector{F}} or a two row Matrix{F} !"));
  end # if
  # NURBS structure represents a surface

  num1 = nurbs.number[1];
  num2 = nurbs.number[2];
  degree = nurbs.order .- 1;

  if typeof(tt)<:Vector{Vector{F}}
    # Evaluate over a [u,v] grid
    # tt[1] represents the u direction
    # tt[2] represents the v direction

    nt1 = length(tt[1]);
    nt2 = length(tt[2]);

    # Evaluate along the v direction
    val = reshape(nurbs.coefs,(4*num1,num2));
    val = bspeval(degree[2],val,nurbs.knots[2],tt[2]);
    val = reshape(val,(4, num1, nt2));

    # Evaluate alon the u direction
    val = permutedims(val,[1,3,2]);
    val = reshape(val,(4*nt2,num1));
    val = bspeval(degree[1],val,nurbs.knots[1],tt[1]);
    val = reshape(val,(4,nt2,nt1));
    val = permutedims(val,[1,3,2]);

    if flag == :cartesian
      w = @view val[4:4,:,:];
      return val[1:3,:,:] ./ repeat(w,outer=[3,1,1]);
    else
      w = val[4:4,:,:];
      p = val[1:3,:,:];
      return p,w;
    end # if
  else
    if typeof(tt)<:Matrix{F} && size(tt,1)!=2
      throw(ArgumentError("Scattered parametric evaluation points for " *
                          "surfaces must be of the type two row Matrix{F} !"));
    end # if
    # Evaluate at scattered points
    # tt[1,:] represents the u direction
    # tt[2,:] represents the v direction

    nt = size(tt,2);

    # Evaluate along the v direction
    val = reshape(nurbs.coefs,(4*num1,num2));
    val = bspeval(degree[2],val,nurbs.knots[2],tt[2,:]);
    val = reshape(val,(4, num1,nt));

    # Evaluate along the u direction
    pnts = zeros(F,(4,nt));
    for v in 1:nt
      pnts[:,v] = bspeval(degree[1],val[:,:,v],nurbs.knots[1],vec([tt[1,v]]));
    end # for v

    if flag == :cartesian
      w = @view pnts[4:4,:];
      return pnts[1:3,:] ./ w;
    else
      w = pnts[4:4,:];
      p = val[1:3,:];
      return p,w;
    end # if
  end # if
else
  if !(typeof(tt)<:Vector{F})
    throw(ArgumentError("Valid parametric evaluation points for curve must " *
                        "be provided as Vector{F} !"));
  end # if
  # NURBS structure represents a curve
  # tt represent a vector of parametric points in the u direction

  val = bspeval(nurbs.order-1,nurbs.coefs,nurbs.knots,tt);

  if flag == :cartesian
    w = @view val[4:4,:];
    return val[1:3,:] ./ w;
  else
    w = val[4:4,:];
    p = val[1:3,:];
    return p,w;
  end # if
end # if
end # nrbeval