diff --git a/Makefile b/Makefile index 2e50708ee..96b0191fb 100644 --- a/Makefile +++ b/Makefile @@ -32,7 +32,6 @@ build_frontend: docs: rm -rf docs/_build/ poetry install --with doc --no-root - # sphinx-apidoc -f -o ./docs/_build/modules ./studio sphinx-autobuild -b html docs docs/_build --port 8001 .PHONY: dockerhub diff --git a/docs/algorithms/caiman.rst b/docs/algorithms/caiman.rst new file mode 100644 index 000000000..6fdf8e82f --- /dev/null +++ b/docs/algorithms/caiman.rst @@ -0,0 +1,23 @@ +CaImAn +====== + +caiman\_cnmf +------------ + +.. automodule:: studio.app.optinist.wrappers.caiman.cnmf + :members: + :undoc-members: + +caiman\_cnmfe +------------- + +.. automodule:: studio.app.optinist.wrappers.caiman.cnmfe + :members: + :undoc-members: + +caiman\_motion\_correction +-------------------------- + +.. automodule:: studio.app.optinist.wrappers.caiman.motion_correction + :members: + :undoc-members: diff --git a/docs/algorithms/index.rst b/docs/algorithms/index.rst new file mode 100644 index 000000000..14ee056ab --- /dev/null +++ b/docs/algorithms/index.rst @@ -0,0 +1,15 @@ +********** +ALGORITHMS +********** + +.. toctree:: + :maxdepth: 2 + + caiman + suite2p + lccd + +.. toctree:: + :maxdepth: 3 + + optinist/index diff --git a/docs/algorithms/lccd.rst b/docs/algorithms/lccd.rst new file mode 100644 index 000000000..b51fb25e7 --- /dev/null +++ b/docs/algorithms/lccd.rst @@ -0,0 +1,9 @@ +Lccd +==== + +lccd\_detection +--------------- + +.. automodule:: studio.app.optinist.wrappers.lccd.lccd_detection + :members: + :undoc-members: diff --git a/docs/algorithms/optinist/basic_neural_analysis.rst b/docs/algorithms/optinist/basic_neural_analysis.rst new file mode 100644 index 000000000..76212618b --- /dev/null +++ b/docs/algorithms/optinist/basic_neural_analysis.rst @@ -0,0 +1,16 @@ +Basic Neural Analysis +===================== + +cell\_grouping +-------------- + +.. automodule:: studio.app.optinist.wrappers.optinist.basic_neural_analysis.cell_grouping + :members: + :undoc-members: + +eta +--- + +.. automodule:: studio.app.optinist.wrappers.optinist.basic_neural_analysis.eta + :members: + :undoc-members: diff --git a/docs/algorithms/optinist/dimension_reduction.rst b/docs/algorithms/optinist/dimension_reduction.rst new file mode 100644 index 000000000..f09b0c873 --- /dev/null +++ b/docs/algorithms/optinist/dimension_reduction.rst @@ -0,0 +1,30 @@ +Dimension Reduction +================== + +cca +--- + +.. automodule:: studio.app.optinist.wrappers.optinist.dimension_reduction.cca + :members: + :undoc-members: + +dpca\_fit +--------- + +.. automodule:: studio.app.optinist.wrappers.optinist.dimension_reduction.dpca_fit + :members: + :undoc-members: + +pca +--- + +.. automodule:: studio.app.optinist.wrappers.optinist.dimension_reduction.pca + :members: + :undoc-members: + +tsne +---- + +.. automodule:: studio.app.optinist.wrappers.optinist.dimension_reduction.tsne + :members: + :undoc-members: diff --git a/docs/algorithms/optinist/index.rst b/docs/algorithms/optinist/index.rst new file mode 100644 index 000000000..39f2fb20b --- /dev/null +++ b/docs/algorithms/optinist/index.rst @@ -0,0 +1,11 @@ +********** +OptiNiSt +********** + +.. toctree:: + :maxdepth: 2 + + basic_neural_analysis + dimension_reduction + neural_decoding + neural_population_analysis diff --git a/docs/algorithms/optinist/neural_decoding.rst b/docs/algorithms/optinist/neural_decoding.rst new file mode 100644 index 000000000..1f30eafa0 --- /dev/null +++ b/docs/algorithms/optinist/neural_decoding.rst @@ -0,0 +1,23 @@ +Neural Decoding +=============== + +glm +--- + +.. automodule:: studio.app.optinist.wrappers.optinist.neural_decoding.glm + :members: + :undoc-members: + +lda +--- + +.. automodule:: studio.app.optinist.wrappers.optinist.neural_decoding.lda + :members: + :undoc-members: + +svm +--- + +.. automodule:: studio.app.optinist.wrappers.optinist.neural_decoding.svm + :members: + :undoc-members: diff --git a/docs/algorithms/optinist/neural_population_analysis.rst b/docs/algorithms/optinist/neural_population_analysis.rst new file mode 100644 index 000000000..edbbdb9d0 --- /dev/null +++ b/docs/algorithms/optinist/neural_population_analysis.rst @@ -0,0 +1,23 @@ +Neural Population Analysis +========================== + +correlation +----------- + +.. automodule:: studio.app.optinist.wrappers.optinist.neural_population_analysis.correlation + :members: + :undoc-members: + +cross\_correlation +------------------ + +.. automodule:: studio.app.optinist.wrappers.optinist.neural_population_analysis.cross_correlation + :members: + :undoc-members: + +granger +------- + +.. automodule:: studio.app.optinist.wrappers.optinist.neural_population_analysis.granger + :members: + :undoc-members: diff --git a/docs/algorithms/suite2p.rst b/docs/algorithms/suite2p.rst new file mode 100644 index 000000000..a19601bd8 --- /dev/null +++ b/docs/algorithms/suite2p.rst @@ -0,0 +1,30 @@ +Suite2p +======= + +suite2p\_file\_convert +---------------------- + +.. automodule:: studio.app.optinist.wrappers.suite2p.file_convert + :members: + :undoc-members: + +suite2p\_registration +--------------------- + +.. automodule:: studio.app.optinist.wrappers.suite2p.registration + :members: + :undoc-members: + +suite2p\_roi +------------ + +.. automodule:: studio.app.optinist.wrappers.suite2p.roi + :members: + :undoc-members: + +suite2p\_spike\_deconv +---------------------- + +.. automodule:: studio.app.optinist.wrappers.suite2p.spike_deconv + :members: + :undoc-members: diff --git a/docs/conf.py b/docs/conf.py index b385f5af0..a1840a6c6 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -5,6 +5,7 @@ # https://www.sphinx-doc.org/en/master/usage/configuration.html import os +import sys from datetime import datetime # -- Path setup -------------------------------------------------------------- @@ -13,7 +14,7 @@ # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # - +sys.path.insert(0, os.path.abspath("../studio")) # -- Project information ----------------------------------------------------- project = "OptiNiSt" @@ -44,11 +45,16 @@ "sphinx.ext.extlinks", "sphinx.ext.autodoc.typehints", "sphinx_copybutton", + "sphinxcontrib.autodoc_pydantic", ] +# autodoc_pydantic config +autodoc_pydantic_model_show_field_summary = False +autodoc_pydantic_settings_show_json_error_strategy = "coerce" +autodoc_pydantic_model_show_json = False + # Tell myst-parser to assign header anchors for h1-h3. myst_heading_anchors = 4 - suppress_warnings = ["myst.header"] # Add any paths that contain templates here, relative to this directory. @@ -79,6 +85,7 @@ html_show_sourcelink = False autosummary_generate = True +napoleon_custom_sections = [("Returns", "params_style")] html_theme_options = { "canonical_url": "", diff --git a/docs/index.rst b/docs/index.rst index a99d40fd2..37885a73a 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -34,6 +34,7 @@ Main Features installation/index tutorials gui/index + algorithms/index host_for_multiuser/index utils/index for_developers/index diff --git a/frontend/src/App.tsx b/frontend/src/App.tsx index 3fdc938ff..8493b88e0 100644 --- a/frontend/src/App.tsx +++ b/frontend/src/App.tsx @@ -33,9 +33,7 @@ const App: FC = () => { dispatch(getModeStandalone()) .unwrap() .catch(() => { - new Promise((resolve) => - setTimeout(resolve, RETRY_WAIT), - ).then(() => { + new Promise((resolve) => setTimeout(resolve, RETRY_WAIT)).then(() => { getMode() }) }) diff --git a/frontend/src/api/nwb/NWB.ts b/frontend/src/api/nwb/NWB.ts index 0a5c32685..de57f2f98 100644 --- a/frontend/src/api/nwb/NWB.ts +++ b/frontend/src/api/nwb/NWB.ts @@ -1,8 +1,8 @@ import { BASE_URL } from "const/API" import axios from "utils/axios" -import { ParamDTO } from "utils/param/ParamType" +import { ParamMap } from "utils/param/ParamType" -export async function getNWBParamsApi(): Promise { +export async function getNWBParamsApi(): Promise { const response = await axios.get(`${BASE_URL}/nwb`) return response.data } diff --git a/frontend/src/api/params/Params.ts b/frontend/src/api/params/Params.ts index 6e2b015fc..7a739ffd0 100644 --- a/frontend/src/api/params/Params.ts +++ b/frontend/src/api/params/Params.ts @@ -1,8 +1,8 @@ import { BASE_URL } from "const/API" import axios from "utils/axios" -import { ParamDTO } from "utils/param/ParamType" +import { ParamMap } from "utils/param/ParamType" -export async function getAlgoParamsApi(algoName: string): Promise { +export async function getAlgoParamsApi(algoName: string): Promise { const response = await axios.get(`${BASE_URL}/params/${algoName}`) return response.data } diff --git a/frontend/src/api/run/Run.ts b/frontend/src/api/run/Run.ts index 11e3ba36e..9e788b782 100644 --- a/frontend/src/api/run/Run.ts +++ b/frontend/src/api/run/Run.ts @@ -7,14 +7,14 @@ import type { } from "store/slice/FlowElement/FlowElementType" import type { FILE_TYPE } from "store/slice/InputNode/InputNodeType" import axios from "utils/axios" -import type { ParamMap } from "utils/param/ParamType" +import type { ParamMapWithoutMeta } from "utils/param/ParamType" export type RunPostData = { name: string nodeDict: NodeDict edgeDict: EdgeDict - nwbParam: ParamMap - snakemakeParam: ParamMap + nwbParam: ParamMapWithoutMeta + snakemakeParam: ParamMapWithoutMeta forceRunList: { nodeId: string; name: string }[] } @@ -39,7 +39,7 @@ export interface InputNodePostData extends InputNodeData { export interface AlgorithmNodePostData extends AlgorithmNodeData { path: string - param: ParamMap + param: ParamMapWithoutMeta } export async function runApi( diff --git a/frontend/src/api/snakemake/Snakemake.ts b/frontend/src/api/snakemake/Snakemake.ts index 01dc28bc4..0f734d550 100644 --- a/frontend/src/api/snakemake/Snakemake.ts +++ b/frontend/src/api/snakemake/Snakemake.ts @@ -1,8 +1,8 @@ import { BASE_URL } from "const/API" import axios from "utils/axios" -import { ParamDTO } from "utils/param/ParamType" +import { ParamMap } from "utils/param/ParamType" -export async function getSnakemakeParamsApi(): Promise { +export async function getSnakemakeParamsApi(): Promise { const response = await axios.get(`${BASE_URL}/snakemake`) return response.data } diff --git a/frontend/src/store/slice/AlgorithmNode/AlgorithmNodeActions.ts b/frontend/src/store/slice/AlgorithmNode/AlgorithmNodeActions.ts index 95c36db66..bb52d1b77 100644 --- a/frontend/src/store/slice/AlgorithmNode/AlgorithmNodeActions.ts +++ b/frontend/src/store/slice/AlgorithmNode/AlgorithmNodeActions.ts @@ -2,10 +2,10 @@ import { createAsyncThunk } from "@reduxjs/toolkit" import { getAlgoParamsApi } from "api/params/Params" import { ALGORITHM_NODE_SLICE_NAME } from "store/slice/AlgorithmNode/AlgorithmNodeType" -import { ParamDTO } from "utils/param/ParamType" +import { ParamMap } from "utils/param/ParamType" export const getAlgoParams = createAsyncThunk< - ParamDTO, + ParamMap, { nodeId: string; algoName: string } >( `${ALGORITHM_NODE_SLICE_NAME}/getAlgoParams`, diff --git a/frontend/src/store/slice/AlgorithmNode/AlgorithmNodeSlice.ts b/frontend/src/store/slice/AlgorithmNode/AlgorithmNodeSlice.ts index 05153a9d9..d858ae985 100644 --- a/frontend/src/store/slice/AlgorithmNode/AlgorithmNodeSlice.ts +++ b/frontend/src/store/slice/AlgorithmNode/AlgorithmNodeSlice.ts @@ -19,7 +19,7 @@ import { importWorkflowConfig, fetchWorkflow, } from "store/slice/Workflow/WorkflowActions" -import { convertToParamMap, getChildParam } from "utils/param/ParamUtils" +import { getChildParam } from "utils/param/ParamUtils" const initialState: AlgorithmNode = {} @@ -50,7 +50,7 @@ export const algorithmNodeSlice = createSlice({ builder .addCase(getAlgoParams.fulfilled, (state, action) => { const { nodeId } = action.meta.arg - state[nodeId].params = convertToParamMap(action.payload) + state[nodeId].params = action.payload }) .addCase(addAlgorithmNode.fulfilled, (state, action) => { const { node, functionPath, name, runAlready } = action.meta.arg @@ -59,7 +59,7 @@ export const algorithmNodeSlice = createSlice({ state[node.id] = { functionPath, name, - params: convertToParamMap(params), + params: params, isUpdated: runAlready ?? false, } } diff --git a/frontend/src/store/slice/FlowElement/FlowElementActions.ts b/frontend/src/store/slice/FlowElement/FlowElementActions.ts index 82834c41c..563b8f153 100644 --- a/frontend/src/store/slice/FlowElement/FlowElementActions.ts +++ b/frontend/src/store/slice/FlowElement/FlowElementActions.ts @@ -8,13 +8,13 @@ import { NodeData, } from "store/slice/FlowElement/FlowElementType" import { FILE_TYPE } from "store/slice/InputNode/InputNodeType" -import { ParamDTO } from "utils/param/ParamType" +import { ParamMap } from "utils/param/ParamType" export type AddingNodeData = Omit, "position"> & Partial, "position">> export const addAlgorithmNode = createAsyncThunk< - ParamDTO, + ParamMap, { node: AddingNodeData functionPath: string diff --git a/frontend/src/store/slice/NWB/NWBAction.ts b/frontend/src/store/slice/NWB/NWBAction.ts index c5041962d..4da8e9d15 100644 --- a/frontend/src/store/slice/NWB/NWBAction.ts +++ b/frontend/src/store/slice/NWB/NWBAction.ts @@ -2,9 +2,9 @@ import { createAsyncThunk } from "@reduxjs/toolkit" import { getNWBParamsApi } from "api/nwb/NWB" import { NWB_SLICE_NAME } from "store/slice/NWB/NWBType" -import { ParamDTO } from "utils/param/ParamType" +import { ParamMap } from "utils/param/ParamType" -export const getNWBParams = createAsyncThunk( +export const getNWBParams = createAsyncThunk( `${NWB_SLICE_NAME}/getNWBParams`, async (_, thunkAPI) => { const { rejectWithValue } = thunkAPI diff --git a/frontend/src/store/slice/NWB/NWBSlice.ts b/frontend/src/store/slice/NWB/NWBSlice.ts index f9d6f414c..125681ae1 100644 --- a/frontend/src/store/slice/NWB/NWBSlice.ts +++ b/frontend/src/store/slice/NWB/NWBSlice.ts @@ -3,7 +3,7 @@ import { createSlice, PayloadAction } from "@reduxjs/toolkit" import { clearFlowElements } from "store/slice/FlowElement/FlowElementSlice" import { getNWBParams } from "store/slice/NWB/NWBAction" import { NWBType, NWB_SLICE_NAME } from "store/slice/NWB/NWBType" -import { convertToParamMap, getChildParam } from "utils/param/ParamUtils" +import { getChildParam } from "utils/param/ParamUtils" const initialState: NWBType = { params: {}, } @@ -29,7 +29,7 @@ export const nwbSlice = createSlice({ extraReducers: (builder) => { builder .addCase(getNWBParams.fulfilled, (state, action) => { - state.params = convertToParamMap(action.payload) + state.params = action.payload }) .addCase(clearFlowElements, () => initialState) }, diff --git a/frontend/src/store/slice/Run/RunSelectors.ts b/frontend/src/store/slice/Run/RunSelectors.ts index faa86c655..1e206a09c 100644 --- a/frontend/src/store/slice/Run/RunSelectors.ts +++ b/frontend/src/store/slice/Run/RunSelectors.ts @@ -32,6 +32,7 @@ import { selectPipelineNodeResultStatus } from "store/slice/Pipeline/PipelineSel import { NODE_RESULT_STATUS } from "store/slice/Pipeline/PipelineType" import { selectSnakemakeParams } from "store/slice/Snakemake/SnakemakeSelectors" import { RootState } from "store/store" +import { removeMetaFromParamMap } from "utils/param/ParamUtils" /** * 前回の結果で、エラーまたはParamに変更があるnodeのリストを返す @@ -57,6 +58,7 @@ const selectNodeDictForRun = (state: RootState): NodeDict => { nodes.forEach((node) => { if (isAlgorithmNodeData(node)) { const param = selectAlgorithmParams(node.id)(state) ?? {} + const paramWithoutMeta = removeMetaFromParamMap(param) const functionPath = selectAlgorithmFunctionPath(node.id)(state) const algorithmNodePostData: Node = { ...node, @@ -65,7 +67,7 @@ const selectNodeDictForRun = (state: RootState): NodeDict => { label: node.data?.label ?? "", type: NODE_TYPE_SET.ALGORITHM, path: functionPath, - param, + param: paramWithoutMeta, }, } nodeDict[node.id] = algorithmNodePostData @@ -113,9 +115,12 @@ export const selectRunPostData = createSelector( nodeDictForRun, forceRunList, ) => { + const nwbParamsWithoutMeta = removeMetaFromParamMap(nwbParams) + const snakemakeParamsWithoutMeta = removeMetaFromParamMap(snakemakeParams) + const runPostData: Omit = { - nwbParam: nwbParams, - snakemakeParam: snakemakeParams, + nwbParam: nwbParamsWithoutMeta, + snakemakeParam: snakemakeParamsWithoutMeta, edgeDict: edgeDictForRun, nodeDict: nodeDictForRun, forceRunList, diff --git a/frontend/src/store/slice/Snakemake/SnakemakeAction.ts b/frontend/src/store/slice/Snakemake/SnakemakeAction.ts index 698bd18e8..c3a0fe241 100644 --- a/frontend/src/store/slice/Snakemake/SnakemakeAction.ts +++ b/frontend/src/store/slice/Snakemake/SnakemakeAction.ts @@ -2,9 +2,9 @@ import { createAsyncThunk } from "@reduxjs/toolkit" import { getSnakemakeParamsApi } from "api/snakemake/Snakemake" import { SNAKEMAKE_SLICE_NAME } from "store/slice/Snakemake/SnakemakeType" -import { ParamDTO } from "utils/param/ParamType" +import { ParamMap } from "utils/param/ParamType" -export const getSnakemakeParams = createAsyncThunk( +export const getSnakemakeParams = createAsyncThunk( `${SNAKEMAKE_SLICE_NAME}/getSnakemakeParams`, async (_, thunkAPI) => { const { rejectWithValue } = thunkAPI diff --git a/frontend/src/store/slice/Snakemake/SnakemakeSlice.ts b/frontend/src/store/slice/Snakemake/SnakemakeSlice.ts index d45603fac..9d9766f0b 100644 --- a/frontend/src/store/slice/Snakemake/SnakemakeSlice.ts +++ b/frontend/src/store/slice/Snakemake/SnakemakeSlice.ts @@ -6,7 +6,7 @@ import { SnakemakeType, SNAKEMAKE_SLICE_NAME, } from "store/slice/Snakemake/SnakemakeType" -import { convertToParamMap, getChildParam } from "utils/param/ParamUtils" +import { getChildParam } from "utils/param/ParamUtils" const initialState: SnakemakeType = { params: {}, @@ -33,7 +33,7 @@ export const SnakemakeSlice = createSlice({ extraReducers: (builder) => { builder .addCase(getSnakemakeParams.fulfilled, (state, action) => { - state.params = convertToParamMap(action.payload) + state.params = action.payload }) .addCase(clearFlowElements, () => initialState) }, diff --git a/frontend/src/store/test/workflow/AlgorithmNode/AlgorithmNode.test.ts b/frontend/src/store/test/workflow/AlgorithmNode/AlgorithmNode.test.ts index 1ee0cc58a..65dbf6724 100644 --- a/frontend/src/store/test/workflow/AlgorithmNode/AlgorithmNode.test.ts +++ b/frontend/src/store/test/workflow/AlgorithmNode/AlgorithmNode.test.ts @@ -22,10 +22,25 @@ describe("AlgorithmNode", () => { const name = "suite2p_roi" const functionPath = "suite2p/suite2p_roi" const algorithmParams = { - param1: 1, - param2: "param2", + param1: { + type: "child", + value: 1, + path: "param1", + }, + param2: { + type: "child", + value: "param2", + path: "param2", + }, param3: { - param4: 4, + type: "parent", + children: { + param4: { + type: "child", + value: 4, + path: "param3/param4", + }, + }, }, } const addAlgorithmNodeAction = { @@ -66,12 +81,12 @@ describe("AlgorithmNode", () => { expect(selectAlgorithmNodeById(nodeId)(targetState).params).toEqual({ param1: { type: "child", - value: algorithmParams.param1, + value: 1, path: "param1", }, param2: { type: "child", - value: algorithmParams.param2, + value: "param2", path: "param2", }, param3: { @@ -79,7 +94,7 @@ describe("AlgorithmNode", () => { children: { param4: { type: "child", - value: algorithmParams.param3.param4, + value: 4, path: "param3/param4", }, }, diff --git a/frontend/src/utils/param/ParamType.ts b/frontend/src/utils/param/ParamType.ts index 0e47bdc19..cf45f502c 100644 --- a/frontend/src/utils/param/ParamType.ts +++ b/frontend/src/utils/param/ParamType.ts @@ -13,10 +13,14 @@ export type ParamParent = { export type ParamChild = { type: "child" + dataType?: string + doc?: string value: unknown path: string } -export type ParamDTO = { - [key: string]: unknown +export type ParamMapWithoutMeta = { + [paramKey: string]: ParamTypeWithoutMeta } +type ParamTypeWithoutMeta = ParamParent | ParamChildWithoutMeta +type ParamChildWithoutMeta = Omit diff --git a/frontend/src/utils/param/ParamUtils.ts b/frontend/src/utils/param/ParamUtils.ts index 43dc58857..cf84fdc56 100644 --- a/frontend/src/utils/param/ParamUtils.ts +++ b/frontend/src/utils/param/ParamUtils.ts @@ -1,9 +1,9 @@ import { - ParamDTO, ParamChild, ParamParent, ParamType, ParamMap, + ParamMapWithoutMeta, } from "utils/param/ParamType" export function getChildParam( @@ -34,71 +34,21 @@ export function isParamParent(param: ParamType): param is ParamParent { return param.type === "parent" } -function isDictObject(value: unknown): value is { [key: string]: unknown } { - return value !== null && typeof value === "object" && !Array.isArray(value) -} - -const PATH_SEPARATOR = "/" - -export function convertToParamMap(dto: ParamDTO, keyList?: string[]): ParamMap { - const ParamMap: ParamMap = {} - Object.entries(dto).forEach(([name, value]) => { - const kList = keyList ?? [] - if (isDictObject(value)) { - kList.push(name) - ParamMap[name] = { - type: "parent", - children: convertToParamMap(value, kList), - } - } else { - ParamMap[name] = { +export function removeMetaFromParamMap(param: ParamMap) { + const result: ParamMapWithoutMeta = {} + for (const [k, v] of Object.entries(param)) { + if (isParamChild(v)) { + result[k] = { type: "child", - value, - path: kList.concat([name]).join(PATH_SEPARATOR), + value: v.value, + path: v.path, + } + } else if (isParamParent(v)) { + result[k] = { + type: "parent", + children: removeMetaFromParamMap(v.children), } } - }) - return ParamMap -} - -export function equalsParamMap(a: ParamMap, b: ParamMap) { - if (a === b) { - return true - } - const aArray = Object.keys(a) - const bArray = Object.keys(b) - return ( - aArray.length === bArray.length && - aArray.every((aKey) => { - const aValue = a[aKey] - const bValue = b[aKey] - return equalsParam(aValue, bValue) - }) - ) -} - -function equalsParam(a: ParamType, b: ParamType): boolean { - if (a === b) { - return true } - if (isParamChild(a) && isParamChild(b)) { - return equalsParamChild(a, b) - } else if (isParamParent(a) && isParamParent(b)) { - const aArray = Object.keys(a) - const bArray = Object.keys(b) - return ( - aArray.length === bArray.length && - aArray.every((aKey) => { - const aValue = a.children[aKey] - const bValue = b.children[aKey] - return equalsParam(aValue, bValue) - }) - ) - } else { - return false - } -} - -function equalsParamChild(a: ParamChild, b: ParamChild) { - return a.path === b.path && a.value === b.value + return result } diff --git a/poetry.lock b/poetry.lock index 78e108c66..1ef38164e 100644 --- a/poetry.lock +++ b/poetry.lock @@ -92,6 +92,27 @@ docs = ["furo", "myst-parser", "sphinx", "sphinx-notfound-page", "sphinxcontrib- tests = ["attrs[tests-no-zope]", "zope-interface"] tests-no-zope = ["cloudpickle", "hypothesis", "mypy (>=1.1.1)", "pympler", "pytest (>=4.3.0)", "pytest-mypy-plugins", "pytest-xdist[psutil]"] +[[package]] +name = "autodoc-pydantic" +version = "1.9.0" +description = "Seamlessly integrate pydantic models in your Sphinx documentation." +optional = false +python-versions = ">=3.7.1,<4.0.0" +files = [ + {file = "autodoc_pydantic-1.9.0-py3-none-any.whl", hash = "sha256:cbf7ec2f27f913629bd38f9944fa6c4a86541c3cadba4a6fa9d2079e500223d8"}, + {file = "autodoc_pydantic-1.9.0.tar.gz", hash = "sha256:0f35f8051abe77b5ae16d8a1084c47a5871435e2ca9060e36c838d063c03cc89"}, +] + +[package.dependencies] +pydantic = ">=1.5,<2.0.0" +Sphinx = ">=3.4" + +[package.extras] +dev = ["coverage (>=7,<8)", "flake8 (>=3,<4)", "pytest (>=7,<8)", "sphinx-copybutton (>=0.4,<0.5)", "sphinx-rtd-theme (>=1.0,<2.0)", "sphinx-tabs (>=3,<4)", "sphinxcontrib-mermaid (>=0.7,<0.8)", "tox (>=3,<4)"] +docs = ["sphinx-copybutton (>=0.4,<0.5)", "sphinx-rtd-theme (>=1.0,<2.0)", "sphinx-tabs (>=3,<4)", "sphinxcontrib-mermaid (>=0.7,<0.8)"] +erdantic = ["erdantic (>=0.5,<0.6)"] +test = ["coverage (>=7,<8)", "pytest (>=7,<8)"] + [[package]] name = "babel" version = "2.13.1" @@ -3070,24 +3091,24 @@ python-versions = ">=3.6" files = [ {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:b42169467c42b692c19cf539c38d4602069d8c1505e97b86387fcf7afb766e1d"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-macosx_13_0_arm64.whl", hash = "sha256:07238db9cbdf8fc1e9de2489a4f68474e70dffcb32232db7c08fa61ca0c7c462"}, + {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-manylinux2014_aarch64.whl", hash = "sha256:d92f81886165cb14d7b067ef37e142256f1c6a90a65cd156b063a43da1708cfd"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl", hash = "sha256:fff3573c2db359f091e1589c3d7c5fc2f86f5bdb6f24252c2d8e539d4e45f412"}, - {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-manylinux_2_24_aarch64.whl", hash = "sha256:aa2267c6a303eb483de8d02db2871afb5c5fc15618d894300b88958f729ad74f"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-musllinux_1_1_i686.whl", hash = "sha256:840f0c7f194986a63d2c2465ca63af8ccbbc90ab1c6001b1978f05119b5e7334"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:024cfe1fc7c7f4e1aff4a81e718109e13409767e4f871443cbff3dba3578203d"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-win32.whl", hash = "sha256:c69212f63169ec1cfc9bb44723bf2917cbbd8f6191a00ef3410f5a7fe300722d"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-win_amd64.whl", hash = "sha256:cabddb8d8ead485e255fe80429f833172b4cadf99274db39abc080e068cbcc31"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:bef08cd86169d9eafb3ccb0a39edb11d8e25f3dae2b28f5c52fd997521133069"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-macosx_13_0_arm64.whl", hash = "sha256:b16420e621d26fdfa949a8b4b47ade8810c56002f5389970db4ddda51dbff248"}, + {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-manylinux2014_aarch64.whl", hash = "sha256:b5edda50e5e9e15e54a6a8a0070302b00c518a9d32accc2346ad6c984aacd279"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl", hash = "sha256:25c515e350e5b739842fc3228d662413ef28f295791af5e5110b543cf0b57d9b"}, - {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-manylinux_2_24_aarch64.whl", hash = "sha256:1707814f0d9791df063f8c19bb51b0d1278b8e9a2353abbb676c2f685dee6afe"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-musllinux_1_1_i686.whl", hash = "sha256:46d378daaac94f454b3a0e3d8d78cafd78a026b1d71443f4966c696b48a6d899"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:09b055c05697b38ecacb7ac50bdab2240bfca1a0c4872b0fd309bb07dc9aa3a9"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-win32.whl", hash = "sha256:53a300ed9cea38cf5a2a9b069058137c2ca1ce658a874b79baceb8f892f915a7"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-win_amd64.whl", hash = "sha256:c2a72e9109ea74e511e29032f3b670835f8a59bbdc9ce692c5b4ed91ccf1eedb"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:ebc06178e8821efc9692ea7544aa5644217358490145629914d8020042c24aa1"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-macosx_13_0_arm64.whl", hash = "sha256:edaef1c1200c4b4cb914583150dcaa3bc30e592e907c01117c08b13a07255ec2"}, + {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-manylinux2014_aarch64.whl", hash = "sha256:7048c338b6c86627afb27faecf418768acb6331fc24cfa56c93e8c9780f815fa"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:d176b57452ab5b7028ac47e7b3cf644bcfdc8cacfecf7e71759f7f51a59e5c92"}, - {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-manylinux_2_24_aarch64.whl", hash = "sha256:1dc67314e7e1086c9fdf2680b7b6c2be1c0d8e3a8279f2e993ca2a7545fecf62"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-musllinux_1_1_i686.whl", hash = "sha256:3213ece08ea033eb159ac52ae052a4899b56ecc124bb80020d9bbceeb50258e9"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:aab7fd643f71d7946f2ee58cc88c9b7bfc97debd71dcc93e03e2d174628e7e2d"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-win32.whl", hash = "sha256:5c365d91c88390c8d0a8545df0b5857172824b1c604e867161e6b3d59a827eaa"}, @@ -3095,7 +3116,7 @@ files = [ {file = "ruamel.yaml.clib-0.2.8-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:a5aa27bad2bb83670b71683aae140a1f52b0857a2deff56ad3f6c13a017a26ed"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:c58ecd827313af6864893e7af0a3bb85fd529f862b6adbefe14643947cfe2942"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-macosx_12_0_arm64.whl", hash = "sha256:f481f16baec5290e45aebdc2a5168ebc6d35189ae6fea7a58787613a25f6e875"}, - {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-manylinux_2_24_aarch64.whl", hash = "sha256:77159f5d5b5c14f7c34073862a6b7d34944075d9f93e681638f6d753606c6ce6"}, + {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-manylinux2014_aarch64.whl", hash = "sha256:3fcc54cb0c8b811ff66082de1680b4b14cf8a81dce0d4fbf665c2265a81e07a1"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:7f67a1ee819dc4562d444bbafb135832b0b909f81cc90f7aa00260968c9ca1b3"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-musllinux_1_1_i686.whl", hash = "sha256:4ecbf9c3e19f9562c7fdd462e8d18dd902a47ca046a2e64dba80699f0b6c09b7"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-musllinux_1_1_x86_64.whl", hash = "sha256:87ea5ff66d8064301a154b3933ae406b0863402a799b16e4a1d24d9fbbcbe0d3"}, @@ -3103,7 +3124,7 @@ files = [ {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-win_amd64.whl", hash = "sha256:3f215c5daf6a9d7bbed4a0a4f760f3113b10e82ff4c5c44bec20a68c8014f675"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:1b617618914cb00bf5c34d4357c37aa15183fa229b24767259657746c9077615"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-macosx_12_0_arm64.whl", hash = "sha256:a6a9ffd280b71ad062eae53ac1659ad86a17f59a0fdc7699fd9be40525153337"}, - {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-manylinux_2_24_aarch64.whl", hash = "sha256:305889baa4043a09e5b76f8e2a51d4ffba44259f6b4c72dec8ca56207d9c6fe1"}, + {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-manylinux2014_aarch64.whl", hash = "sha256:665f58bfd29b167039f714c6998178d27ccd83984084c286110ef26b230f259f"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:700e4ebb569e59e16a976857c8798aee258dceac7c7d6b50cab63e080058df91"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-musllinux_1_1_i686.whl", hash = "sha256:e2b4c44b60eadec492926a7270abb100ef9f72798e18743939bdbf037aab8c28"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:e79e5db08739731b0ce4850bed599235d601701d5694c36570a99a0c5ca41a9d"}, @@ -3111,7 +3132,7 @@ files = [ {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-win_amd64.whl", hash = "sha256:56f4252222c067b4ce51ae12cbac231bce32aee1d33fbfc9d17e5b8d6966c312"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:03d1162b6d1df1caa3a4bd27aa51ce17c9afc2046c31b0ad60a0a96ec22f8001"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-macosx_12_0_arm64.whl", hash = "sha256:bba64af9fa9cebe325a62fa398760f5c7206b215201b0ec825005f1b18b9bccf"}, - {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-manylinux_2_24_aarch64.whl", hash = "sha256:a1a45e0bb052edf6a1d3a93baef85319733a888363938e1fc9924cb00c8df24c"}, + {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-manylinux2014_aarch64.whl", hash = "sha256:9eb5dee2772b0f704ca2e45b1713e4e5198c18f515b52743576d196348f374d3"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:da09ad1c359a728e112d60116f626cc9f29730ff3e0e7db72b9a2dbc2e4beed5"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-musllinux_1_1_i686.whl", hash = "sha256:184565012b60405d93838167f425713180b949e9d8dd0bbc7b49f074407c5a8b"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:a75879bacf2c987c003368cf14bed0ffe99e8e85acfa6c0bfffc21a090f16880"}, @@ -3980,4 +4001,4 @@ testing = ["big-O", "jaraco.functools", "jaraco.itertools", "more-itertools", "p [metadata] lock-version = "2.0" python-versions = "^3.8.1" -content-hash = "24a03c16411f6e968310732e134badc23d1693f98e6c5e407e722637aa64f804" +content-hash = "4f0820d4cc75ce68d6213076dd2ea30580995474ecd22bd41c128845c98d834e" diff --git a/pyproject.toml b/pyproject.toml index b551ecd36..695dcf2b3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,8 +28,6 @@ include = [ "frontend/build/static/css/*", "frontend/build/static/js/*", "frontend/build/static/media/*", - "studio/app/optinist/core/nwb/nwb.yaml", - "studio/app/common/core/snakemake/snakemake.yaml", "studio/app/*/wrappers/**/*.yaml", "studio/app/Snakefile", "studio/config/*.yaml", @@ -91,6 +89,7 @@ sphinx-autodoc-typehints = "*" sphinx-copybutton = "0.5.2" sphinx-autobuild = "2021.3.14" myst-parser = "*" +autodoc-pydantic = "1.9.0" [tool.poetry.group.test] optional = true diff --git a/studio/app/common/core/algo.py b/studio/app/common/core/algo.py new file mode 100644 index 000000000..8234efc63 --- /dev/null +++ b/studio/app/common/core/algo.py @@ -0,0 +1,16 @@ +from abc import ABC, abstractmethod + + +class AlgoTemplate(ABC): + def set_output_dir(self, output_dir: str): + self.output_dir = output_dir + self.function_id = output_dir.split("/")[-1] + return self + + def set_nwb_params(self, nwb_params): + self.nwb_params = nwb_params + return self + + @abstractmethod + def run(self, params): + pass diff --git a/studio/app/common/core/rules/runner.py b/studio/app/common/core/rules/runner.py index 408a95845..4f0742fa7 100644 --- a/studio/app/common/core/rules/runner.py +++ b/studio/app/common/core/rules/runner.py @@ -123,11 +123,14 @@ def save_all_nwb(cls, save_path, all_nwbfile): @classmethod def execute_function(cls, path, params, nwb_params, output_dir, input_info): wrapper = cls.dict2leaf(wrapper_dict, path.split("/")) - func = copy.deepcopy(wrapper["function"]) - output_info = func( - params=params, nwbfile=nwb_params, output_dir=output_dir, **input_info + algo = copy.deepcopy(wrapper["function"]) + output_info = ( + algo() + .set_output_dir(output_dir) + .set_nwb_params(nwb_params) + .run(params=params, **input_info) ) - del func + del algo gc.collect() return output_info diff --git a/studio/app/common/core/snakemake/snakemake.yaml b/studio/app/common/core/snakemake/snakemake.yaml deleted file mode 100644 index 3de99490b..000000000 --- a/studio/app/common/core/snakemake/snakemake.yaml +++ /dev/null @@ -1,5 +0,0 @@ -use_conda: True -cores: 2 -forceall: False -forcetargets: True -lock: False diff --git a/studio/app/common/core/snakemake/snakemake_rule.py b/studio/app/common/core/snakemake/snakemake_rule.py index 79b38c711..d87c4ad7b 100644 --- a/studio/app/common/core/snakemake/snakemake_rule.py +++ b/studio/app/common/core/snakemake/snakemake_rule.py @@ -3,8 +3,8 @@ from studio.app.common.core.snakemake.smk import Rule from studio.app.common.core.snakemake.smk_builder import RuleBuilder from studio.app.common.core.utils.filepath_creater import get_pickle_file +from studio.app.common.core.utils.param_utils import get_typecheck_params from studio.app.common.core.workflow.workflow import Edge, Node, NodeType -from studio.app.common.core.workflow.workflow_params import get_typecheck_params class SmkRule: diff --git a/studio/app/common/core/utils/filepath_finder.py b/studio/app/common/core/utils/filepath_finder.py index 31c19ab4d..f17bb37d8 100644 --- a/studio/app/common/core/utils/filepath_finder.py +++ b/studio/app/common/core/utils/filepath_finder.py @@ -3,7 +3,7 @@ from typing import Optional from studio.app.common.core.utils.filepath_creater import join_filepath -from studio.app.dir_path import CORE_PARAM_PATH, DIRPATH +from studio.app.dir_path import DIRPATH def find_filepath(name, category) -> Optional[str]: @@ -18,12 +18,5 @@ def find_filepath(name, category) -> Optional[str]: return filepaths[0] if len(filepaths) > 0 else None -def find_param_filepath(name: str): - if name in CORE_PARAM_PATH.__members__: - return CORE_PARAM_PATH[name].value - else: - return find_filepath(name, "params") - - def find_condaenv_filepath(name: str): return find_filepath(name, "conda") diff --git a/studio/app/common/core/utils/param_utils.py b/studio/app/common/core/utils/param_utils.py new file mode 100644 index 000000000..a254b8d21 --- /dev/null +++ b/studio/app/common/core/utils/param_utils.py @@ -0,0 +1,83 @@ +import inspect +from typing import Any, Dict, _GenericAlias + +from pydantic import BaseModel + +from studio.app.common.core.algo import AlgoTemplate +from studio.app.common.schemas.params import SnakemakeParams +from studio.app.optinist.schemas.nwb import NWBParams +from studio.app.wrappers import wrapper_dict + + +def find_algo_by_name( + name: str, + wrapper_dict: Dict[str, Any] = wrapper_dict, +) -> AlgoTemplate: + for k, v in wrapper_dict.items(): + if isinstance(v, dict) and "function" not in v: + algo = find_algo_by_name(name, v) + if algo is not None: + return algo + elif k == name: + return v["function"] + return None + + +def get_default_params(name: str): + if name == "nwb": + return NWBParams + elif name == "snakemake": + return SnakemakeParams + else: + algo = find_algo_by_name(name) + if algo is None: + return None + else: + sig = inspect.signature(algo.run) + return sig.parameters.get("params").annotation + + +def get_param_map(params, path=None): + default_params = {} + + for k, v in params.__fields__.items(): + k = k if v.field_info.alias is None else v.field_info.alias + if issubclass(type(v.default), BaseModel): + default_params[k] = { + "type": "parent", + "children": get_param_map( + params=v.default, + path=k if path is None else f"{path}/{k}", + ), + } + else: + default_params[k] = { + "type": "child", + "dataType": type(v.default).__name__ + if isinstance(v.annotation, _GenericAlias) + else v.annotation.__name__, + "value": v.default, + "doc": v.field_info.description, + "path": k if path is None else f"{path}/{k}", + } + + return default_params + + +def get_typecheck_params(message_params, name): + default_params = get_default_params(name) + if message_params != {} and message_params is not None: + return default_params(**nest2dict(message_params)).dict(by_alias=True) + else: + return default_params().dict(by_alias=True) + + +def nest2dict(value): + nwb_dict = {} + for _k, _v in value.items(): + if _v["type"] == "child": + nwb_dict[_k] = _v["value"] + elif _v["type"] == "parent": + nwb_dict[_k] = nest2dict(_v["children"]) + + return nwb_dict diff --git a/studio/app/common/core/workflow/workflow_params.py b/studio/app/common/core/workflow/workflow_params.py deleted file mode 100644 index 88e67627d..000000000 --- a/studio/app/common/core/workflow/workflow_params.py +++ /dev/null @@ -1,38 +0,0 @@ -from studio.app.common.core.utils.config_handler import ConfigReader -from studio.app.common.core.utils.filepath_finder import find_param_filepath - - -def get_typecheck_params(message_params, name): - default_params = ConfigReader.read(find_param_filepath(name)) - if message_params != {} and message_params is not None: - return check_types(nest2dict(message_params), default_params) - return default_params - - -def check_types(params, default_params): - for key in params.keys(): - if isinstance(params[key], dict): - params[key] = check_types(params[key], default_params[key]) - else: - if not isinstance(type(params[key]), type(default_params[key])): - data_type = type(default_params[key]) - p = params[key] - if isinstance(data_type, str): - params[key] = str(p) - elif isinstance(data_type, float): - params[key] = float(p) - elif isinstance(data_type, int): - params[key] = int(p) - - return params - - -def nest2dict(value): - nwb_dict = {} - for _k, _v in value.items(): - if _v["type"] == "child": - nwb_dict[_k] = _v["value"] - elif _v["type"] == "parent": - nwb_dict[_k] = nest2dict(_v["children"]) - - return nwb_dict diff --git a/studio/app/common/core/workflow/workflow_result.py b/studio/app/common/core/workflow/workflow_result.py index dd649c0b5..5cd346ed6 100644 --- a/studio/app/common/core/workflow/workflow_result.py +++ b/studio/app/common/core/workflow/workflow_result.py @@ -88,31 +88,37 @@ def cancel(self): """ The algorithm function of this workflow is being executed at the line: https://github.com/snakemake/snakemake/blob/27b224ed12448df8aebc7d1ff8f25e3bf7622232/snakemake/shell.py#L258 - ``` - proc = sp.Popen( - cmd, - bufsize=-1, - shell=use_shell, - stdout=stdout, - universal_newlines=iterable or read or None, - close_fds=close_fds, - **cls._process_args, - env=envvars, - ) - ``` + + :: + + proc = sp.Popen( + cmd, + bufsize=-1, + shell=use_shell, + stdout=stdout, + universal_newlines=iterable or read or None, + close_fds=close_fds, + **cls._process_args, + env=envvars, + ) + The `cmd` argument has the following format: - ``` - source ~/miniconda3/bin/activate - '~/Documents/optinistfs/.snakemake/conda/491889952d2f07f3876bb801eea626e9_'; - set -euo pipefail; - python ~/Documents/optinistfs/.snakemake/scripts/tmp03froqxo.func.py - ``` + + :: + + source ~/miniconda3/bin/activate + '~/Documents/optinistfs/.snakemake/conda/491889952d2f07f3876bb801eea626e9_'; + set -euo pipefail; + python ~/Documents/optinistfs/.snakemake/scripts/tmp03froqxo.func.py + Interrupt the conda activate at the beginning of the process is impossible because it is only called when each algorithm function executes. This workflow is cancelled by killing process via PID of algorithm function saved in pid.json file + Raises: HTTPException: if pid_filepath or last_script_file does not exist + """ if not os.path.exists(self.pid_filepath): raise HTTPException(status_code=404) diff --git a/studio/app/common/core/workflow/workflow_runner.py b/studio/app/common/core/workflow/workflow_runner.py index f72a2d877..e2f91fa41 100644 --- a/studio/app/common/core/workflow/workflow_runner.py +++ b/studio/app/common/core/workflow/workflow_runner.py @@ -10,8 +10,8 @@ from studio.app.common.core.snakemake.snakemake_reader import SmkParamReader from studio.app.common.core.snakemake.snakemake_rule import SmkRule from studio.app.common.core.snakemake.snakemake_writer import SmkConfigWriter +from studio.app.common.core.utils.param_utils import get_typecheck_params from studio.app.common.core.workflow.workflow import NodeType, RunItem -from studio.app.common.core.workflow.workflow_params import get_typecheck_params from studio.app.common.core.workflow.workflow_reader import WorkflowConfigReader from studio.app.common.core.workflow.workflow_writer import WorkflowConfigWriter diff --git a/studio/app/common/routers/algolist.py b/studio/app/common/routers/algolist.py index b5526c173..9c78d7ebb 100644 --- a/studio/app/common/routers/algolist.py +++ b/studio/app/common/routers/algolist.py @@ -21,7 +21,7 @@ def get_nest_dict(cls, parent_value, parent_key: str) -> Dict[str, Algo]: value, cls._parent_key(parent_key, key) ) else: - sig = inspect.signature(value["function"]) + sig = inspect.signature(value["function"].run) returns_list = None if sig.return_annotation is not inspect._empty: returns_list = cls._return_list(sig.return_annotation.items()) @@ -29,7 +29,6 @@ def get_nest_dict(cls, parent_value, parent_key: str) -> Dict[str, Algo]: algo_dict[key] = Algo( args=cls._args_list(sig.parameters.values()), returns=returns_list, - parameter=value["parameter"] if "parameter" in value else None, path=cls._parent_key(parent_key, key), ) @@ -63,23 +62,25 @@ def _parent_key(cls, parent_key: str, key: str) -> str: async def get_algolist() -> Dict[str, Algo]: """_summary_ - Returns: - { - 'caiman': { - 'children': { - 'caiman_mc' : { - 'args': ['images', 'timeseries'], - 'return': ['images'], - 'path': 'caiman/caiman_mc' - }, - 'caiman_cnmf': { - 'args': ['images', 'timeseries'], - 'return': ['images'], - 'path': 'caiman/caiman_mc' + Returns:: + + { + 'caiman': { + 'children': { + 'caiman_mc' : { + 'args': ['images', 'timeseries'], + 'return': ['images'], + 'path': 'caiman/caiman_mc' + }, + 'caiman_cnmf': { + 'args': ['images', 'timeseries'], + 'return': ['images'], + 'path': 'caiman/caiman_mc' + } } } } - } + """ return NestDictGetter.get_nest_dict(wrapper_dict, "") diff --git a/studio/app/common/routers/params.py b/studio/app/common/routers/params.py index ec1eeda78..715cc4101 100644 --- a/studio/app/common/routers/params.py +++ b/studio/app/common/routers/params.py @@ -1,22 +1,22 @@ -from typing import Any, Dict +from typing import Dict, Union -from fastapi import APIRouter +from fastapi import APIRouter, HTTPException -from studio.app.common.core.utils.config_handler import ConfigReader -from studio.app.common.core.utils.filepath_finder import find_param_filepath -from studio.app.common.schemas.params import SnakemakeParams +from studio.app.common.core.utils.param_utils import get_default_params, get_param_map +from studio.app.common.schemas.params import ParamChild, ParamParent, SnakemakeParams router = APIRouter(tags=["params"]) -@router.get("/params/{name}", response_model=Dict[str, Any]) -async def get_params(name: str): - filepath = find_param_filepath(name) - config = ConfigReader.read(filepath) - return config +@router.get("/params/{name}", response_model=Dict[str, Union[ParamChild, ParamParent]]) +async def get_params(name: str) -> dict: + default_params = get_default_params(name) + if default_params is None: + raise HTTPException(status_code=404, detail=f"Algo with name {name} not found") + return get_param_map(default_params) -@router.get("/snakemake", response_model=SnakemakeParams) + +@router.get("/snakemake", response_model=Dict[str, Union[ParamChild, ParamParent]]) async def get_snakemake_params(): - filepath = find_param_filepath("snakemake") - return ConfigReader.read(filepath) + return get_param_map(SnakemakeParams) diff --git a/studio/app/common/schemas/algolist.py b/studio/app/common/schemas/algolist.py index e97a0f1aa..d01e14d2b 100644 --- a/studio/app/common/schemas/algolist.py +++ b/studio/app/common/schemas/algolist.py @@ -21,7 +21,6 @@ class Return: class Algo: args: List[Arg] returns: List[Return] - parameter: str = None path: str = None @@ -36,7 +35,6 @@ class AlgoList(BaseModel): "caiman_mc": { "args": [{"name": "image", "type": "ImageData", "isNone": False}], "returns": [{"name": "mc_images", "type": "ImageData"}], - "parameter": None, "path": "caiman/caiman_mc", } } diff --git a/studio/app/common/schemas/params.py b/studio/app/common/schemas/params.py index d66736e39..8dea6ee07 100644 --- a/studio/app/common/schemas/params.py +++ b/studio/app/common/schemas/params.py @@ -1,9 +1,25 @@ +from typing import Any, Dict, Optional, Union + from pydantic import BaseModel +from pydantic.dataclasses import Field + + +class ParamChild(BaseModel): + type: str = "child" + dataType: str + value: Any + path: str + doc: Optional[str] = None + + +class ParamParent(BaseModel): + type: str = "parent" + children: Dict[str, Union[ParamChild, "ParamParent"]] class SnakemakeParams(BaseModel): - use_conda: bool - cores: int - forceall: bool - forcetargets: bool - lock: bool + use_conda: bool = Field(True) + cores: int = Field(2) + forceall: bool = Field(False) + forcetargets: bool = Field(True) + lock: bool = Field(False) diff --git a/studio/app/const.py b/studio/app/const.py index 0192415f1..c24cbc559 100644 --- a/studio/app/const.py +++ b/studio/app/const.py @@ -13,6 +13,6 @@ class FILETYPE: ACCEPT_CSV_EXT = [".csv"] ACCEPT_HDF5_EXT = [".hdf5", ".nwb", ".HDF5", ".NWB"] -NOT_DISPLAY_ARGS_LIST = ["params", "output_dir", "nwbfile", "kwargs"] +NOT_DISPLAY_ARGS_LIST = ["self", "params", "kwargs"] DATE_FORMAT = "%Y-%m-%d %H:%M:%S" diff --git a/studio/app/dir_path.py b/studio/app/dir_path.py index edaac3c2b..2a858d75b 100644 --- a/studio/app/dir_path.py +++ b/studio/app/dir_path.py @@ -1,5 +1,4 @@ import os -from enum import Enum _DEFAULT_DIR = "/tmp/studio" _ENV_DIR = os.environ.get("OPTINIST_DIR") @@ -35,8 +34,3 @@ class DIRPATH: FIREBASE_PRIVATE_PATH = f"{CONFIG_DIR}/auth/firebase_private.json" FIREBASE_CONFIG_PATH = f"{CONFIG_DIR}/auth/firebase_config.json" - - -class CORE_PARAM_PATH(Enum): - nwb = f"{DIRPATH.APP_DIR}/optinist/core/nwb/nwb.yaml" - snakemake = f"{DIRPATH.APP_DIR}/common/core/snakemake/snakemake.yaml" diff --git a/studio/app/optinist/core/nwb/nwb.yaml b/studio/app/optinist/core/nwb/nwb.yaml deleted file mode 100644 index 9cb41ed90..000000000 --- a/studio/app/optinist/core/nwb/nwb.yaml +++ /dev/null @@ -1,25 +0,0 @@ -session_description: 'optinist' -identifier: 'optinist' -experiment_description: 'None' -device: - name: 'Microscope device' - description: 'Microscope Information' - manufacturer: 'Microscope Manufacture' -optical_channel: - name: 'OpticalChannel' - description: 'optical channel' - emission_lambda: 500.0 -imaging_plane: - name: 'ImagingPlane' - description: 'standard' - imaging_rate: 30.0 - excitation_lambda: 900.0 - indicator: 'GCaMP' - location: 'V1' -image_series: - starting_time: 0 - starting_frame: [0,] -ophys: - plane_segmentation: - name: 'PlaneSegmentation' - description: '' diff --git a/studio/app/optinist/routers/nwb.py b/studio/app/optinist/routers/nwb.py index ff172be5b..41cf638e3 100644 --- a/studio/app/optinist/routers/nwb.py +++ b/studio/app/optinist/routers/nwb.py @@ -1,24 +1,26 @@ from glob import glob +from typing import Dict, Union from fastapi import APIRouter, Depends from fastapi.responses import FileResponse -from studio.app.common.core.utils.config_handler import ConfigReader from studio.app.common.core.utils.filepath_creater import join_filepath -from studio.app.common.core.utils.filepath_finder import find_param_filepath +from studio.app.common.core.utils.param_utils import get_param_map from studio.app.common.core.workspace.workspace_dependencies import ( is_workspace_available, ) +from studio.app.common.schemas.params import ParamChild, ParamParent from studio.app.dir_path import DIRPATH from studio.app.optinist.schemas.nwb import NWBParams router = APIRouter() -@router.get("/nwb", response_model=NWBParams, tags=["params"]) +@router.get( + "/nwb", tags=["params"], response_model=Dict[str, Union[ParamChild, ParamParent]] +) async def get_nwb_params(): - filepath = find_param_filepath("nwb") - return ConfigReader.read(filepath) + return get_param_map(NWBParams) @router.get( diff --git a/studio/app/optinist/schemas/nwb.py b/studio/app/optinist/schemas/nwb.py index 2f51e0786..70bc0868c 100644 --- a/studio/app/optinist/schemas/nwb.py +++ b/studio/app/optinist/schemas/nwb.py @@ -1,14 +1,54 @@ -from typing import Any, Dict, Optional, Union +from typing import Optional from pydantic import BaseModel +from pydantic.dataclasses import Field + + +class Device(BaseModel): + name: str = Field("Microscope device") + description: str = Field("Microscope Information") + manufacturer: str = Field("Microscope Manufacture") + + +class OpticalChannel(BaseModel): + name: str = Field("OpticalChannel") + description: str = Field("optical channel") + emission_lambda: float = Field(500.0) + + +class ImagingPlane(BaseModel): + name: str = Field("ImagingPlane") + description: str = Field("standard") + imaging_rate: float = Field(30.0) + excitation_lambda: float = Field(900.0) + indicator: str = Field("GCaMP") + location: str = Field("V1") + + +class ImageSeries(BaseModel): + starting_time: int = Field(0) + starting_frame: list = Field( + [ + 0, + ] + ) + + +class PlaneSegmentation(BaseModel): + name: str = Field("PlaneSegmentation") + description: str = Field("") + + +class Ophys(BaseModel): + plane_segmentation: PlaneSegmentation = Field(PlaneSegmentation()) class NWBParams(BaseModel): - session_description: str = "optinist" - identifier: str = "optinist" - experiment_description: Optional[str] = None - device: Union[Dict, Any] - optical_channel: Union[Dict, Any] - imaging_plane: Union[Dict, Any] - image_series: Union[Dict, Any] - ophys: Union[Dict, Any] + session_description: str = Field("optinist") + identifier: str = Field("optinist") + experiment_description: Optional[str] = Field(None) + device: Device = Field(Device()) + optical_channel: OpticalChannel = Field(OpticalChannel()) + imaging_plane: ImagingPlane = Field(ImagingPlane()) + image_series: ImageSeries = Field(ImageSeries()) + ophys: Ophys = Field(Ophys()) diff --git a/studio/app/optinist/wrappers/caiman/__init__.py b/studio/app/optinist/wrappers/caiman/__init__.py index e0a3d5cc1..876ec9577 100644 --- a/studio/app/optinist/wrappers/caiman/__init__.py +++ b/studio/app/optinist/wrappers/caiman/__init__.py @@ -1,19 +1,19 @@ -from studio.app.optinist.wrappers.caiman.cnmf import caiman_cnmf -from studio.app.optinist.wrappers.caiman.cnmfe import caiman_cnmfe -from studio.app.optinist.wrappers.caiman.motion_correction import caiman_mc +from studio.app.optinist.wrappers.caiman.cnmf import CaimanCnmf +from studio.app.optinist.wrappers.caiman.cnmfe import CaimanCnmfE +from studio.app.optinist.wrappers.caiman.motion_correction import CaimanMc caiman_wrapper_dict = { "caiman": { "caiman_mc": { - "function": caiman_mc, + "function": CaimanMc, "conda_name": "caiman", }, "caiman_cnmf": { - "function": caiman_cnmf, + "function": CaimanCnmf, "conda_name": "caiman", }, "caiman_cnmfe": { - "function": caiman_cnmfe, + "function": CaimanCnmfE, "conda_name": "caiman", }, } diff --git a/studio/app/optinist/wrappers/caiman/cnmf.py b/studio/app/optinist/wrappers/caiman/cnmf.py index 13e20aa7f..0f31af9ee 100644 --- a/studio/app/optinist/wrappers/caiman/cnmf.py +++ b/studio/app/optinist/wrappers/caiman/cnmf.py @@ -1,266 +1,366 @@ import gc +from typing import List, Optional import numpy as np +from pydantic import BaseModel +from pydantic.dataclasses import Field +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.core.utils.filepath_creater import join_filepath from studio.app.common.dataclass import ImageData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import CaimanCnmfData, FluoData, IscellData, RoiData -def get_roi(A, thr, thr_method, swap_dim, dims): - from scipy.ndimage import binary_fill_holes - from skimage.measure import find_contours - - d, nr = np.shape(A) - - # for each patches - ims = [] - coordinates = [] - for i in range(nr): - pars = dict() - # we compute the cumulative sum of the energy of the Ath component - # that has been ordered from least to highest - patch_data = A.data[A.indptr[i] : A.indptr[i + 1]] - indx = np.argsort(patch_data)[::-1] - - if thr_method == "nrg": - cumEn = np.cumsum(patch_data[indx] ** 2) - if len(cumEn) == 0: - pars = dict( - coordinates=np.array([]), - CoM=np.array([np.NaN, np.NaN]), - neuron_id=i + 1, - ) - coordinates.append(pars) - continue - else: - # we work with normalized values - cumEn /= cumEn[-1] - Bvec = np.ones(d) - # we put it in a similar matrix - Bvec[A.indices[A.indptr[i] : A.indptr[i + 1]][indx]] = cumEn - else: - Bvec = np.zeros(d) - Bvec[A.indices[A.indptr[i] : A.indptr[i + 1]]] = ( - patch_data / patch_data.max() - ) +class CaimanCnmfInitParams(BaseModel): + # Ain: Optional[List] = Field( + # None, description="possibility to seed with predetermined binary masks" + # ) + do_refit: bool = Field( + False, + description=( + "Whether to re-run seeded CNMF on accepted patches " + "to refine and perform deconvolution" + ), + ) + K: Optional[int] = Field( + 4, description="upper bound on number of components per patch, in general None" + ) + gSig: List[int] = Field( + [4, 4], + description=( + "gaussian width of a 2D gaussian kernel, which approximates a neuron" + ), + ) + ssub: int = Field( + 1, + description=( + "downsampling factor in space for initialization, " + "increase if you have memory problems. " + "you can pass them here as boolean vectors" + ), + ) + tsub: int = Field( + 1, + description=( + "downsampling factor in time for initialization, " + "increase if you have memory problems." + ), + ) + nb: int = Field( + 2, + description=""" + number of background components (rank) if positive, + else exact ring model with following settings. - if swap_dim: - Bmat = np.reshape(Bvec, dims, order="C") - else: - Bmat = np.reshape(Bvec, dims, order="F") - - r_mask = np.zeros_like(Bmat, dtype="bool") - contour = find_contours(Bmat, thr) - for c in contour: - r_mask[np.round(c[:, 0]).astype("int"), np.round(c[:, 1]).astype("int")] = 1 - - # Fill in the hole created by the contour boundary - r_mask = binary_fill_holes(r_mask) - ims.append(r_mask + (i * r_mask)) - - return ims - - -def caiman_cnmf( - images: ImageData, output_dir: str, params: dict = None, **kwargs -) -> dict(fluorescence=FluoData, iscell=IscellData): - import scipy - from caiman import local_correlations, stop_server - from caiman.cluster import setup_cluster - from caiman.mmapping import prepare_shape - from caiman.paths import memmap_frames_filename - from caiman.source_extraction.cnmf import cnmf - from caiman.source_extraction.cnmf.params import CNMFParams - - function_id = output_dir.split("/")[-1] - print("start caiman_cnmf:", function_id) - - # flatten params segments. - params_flatten = {} - for params_segment in params.values(): - params_flatten.update(params_segment) - params = params_flatten - - Ain = params.pop("Ain", None) - do_refit = params.pop("do_refit", None) - thr = params.pop("thr", None) - - file_path = images.path - if isinstance(file_path, list): - file_path = file_path[0] - - images = images.data - - # np.arrayをmmapへ変換 - order = "C" - dims = images.shape[1:] - T = images.shape[0] - shape_mov = (np.prod(dims), T) - - dir_path = join_filepath(file_path.split("/")[:-1]) - basename = file_path.split("/")[-1] - fname_tot = memmap_frames_filename(basename, dims, T, order) - - mmap_images = np.memmap( - join_filepath([dir_path, fname_tot]), - mode="w+", - dtype=np.float32, - shape=prepare_shape(shape_mov), - order=order, + gnb= 0: Return background as b and W + + gnb=-1: Return full rank background B + + gnb<-1: Don't return background + """, ) + method_init: str = Field("greedy_roi") - mmap_images = np.reshape(mmap_images.T, [T] + list(dims), order="F") - mmap_images[:] = images[:] - del images - gc.collect() +class CaimanCnmfPreprocessParams(BaseModel): + p: int = Field(1, description="order of the autoregressive system") - nwbfile = kwargs.get("nwbfile", {}) - fr = nwbfile.get("imaging_plane", {}).get("imaging_rate", 30) - if params is None: - ops = CNMFParams() - else: - ops = CNMFParams(params_dict={**params, "fr": fr}) +class CaimanCnmfPatchParams(BaseModel): + rf: Optional[int] = Field( + None, + description=( + "half-size of the patches in pixels. e.g., if rf=40, patches are 80x80" + ), + ) + stride: int = Field( + 6, + description=( + "amount of overlap between the patches in pixels, " + "(keep it at least large as gSiz, i.e 4 times the neuron size gSig)" + ), + ) - if "dview" in locals(): - stop_server(dview=dview) # noqa: F821 - c, dview, n_processes = setup_cluster( - backend="local", n_processes=None, single_thread=True +class CaimanCnmfMergeParams(BaseModel): + thr: float = Field(0.9) + merge_thr: float = Field( + 0.85, description="merging threshold, max correlation allowed" ) - cnm = cnmf.CNMF(n_processes=n_processes, dview=dview, Ain=Ain, params=ops) - cnm = cnm.fit(mmap_images) - if do_refit: - cnm = cnm.refit(mmap_images, dview=dview) +class CaimanCnmfParams(BaseModel): + init_params: CaimanCnmfInitParams = Field(CaimanCnmfInitParams()) + preprocess_params: CaimanCnmfPreprocessParams = Field(CaimanCnmfPreprocessParams()) + patch_params: CaimanCnmfPatchParams = Field(CaimanCnmfPatchParams()) + merge_params: CaimanCnmfMergeParams = Field(CaimanCnmfMergeParams()) + + +class CaimanCnmf(AlgoTemplate): + def run( + self, params: CaimanCnmfParams, images: ImageData + ) -> dict(fluorescence=FluoData, iscell=IscellData): + import scipy + from caiman import local_correlations, stop_server + from caiman.cluster import setup_cluster + from caiman.mmapping import prepare_shape + from caiman.paths import memmap_frames_filename + from caiman.source_extraction.cnmf import cnmf + from caiman.source_extraction.cnmf.params import CNMFParams + + print("start caiman_cnmf:", self.function_id) + + # flatten params segments. + params_flatten = {} + for params_segment in params.values(): + params_flatten.update(params_segment) + params = params_flatten + + Ain = params.pop("Ain", None) + do_refit = params.pop("do_refit", None) + thr = params.pop("thr", None) + + file_path = images.path + if isinstance(file_path, list): + file_path = file_path[0] + + images = images.data + + # np.arrayをmmapへ変換 + order = "C" + dims = images.shape[1:] + T = images.shape[0] + shape_mov = (np.prod(dims), T) + + dir_path = join_filepath(file_path.split("/")[:-1]) + basename = file_path.split("/")[-1] + fname_tot = memmap_frames_filename(basename, dims, T, order) + + mmap_images = np.memmap( + join_filepath([dir_path, fname_tot]), + mode="w+", + dtype=np.float32, + shape=prepare_shape(shape_mov), + order=order, + ) - stop_server(dview=dview) + mmap_images = np.reshape(mmap_images.T, [T] + list(dims), order="F") + mmap_images[:] = images[:] - # contours plot - Cn = local_correlations(mmap_images.transpose(1, 2, 0)) - Cn[np.isnan(Cn)] = 0 + del images + gc.collect() - thr_method = "nrg" - swap_dim = False + fr = self.nwb_params.get("imaging_plane", {}).get("imaging_rate", 30) - iscell = np.concatenate( - [ - np.ones(cnm.estimates.A.shape[-1]), - np.zeros(cnm.estimates.b.shape[-1] if cnm.estimates.b is not None else 0), - ] - ).astype(bool) + if params is None: + ops = CNMFParams() + else: + ops = CNMFParams(params_dict={**params, "fr": fr}) - ims = get_roi(cnm.estimates.A, thr, thr_method, swap_dim, dims) - ims = np.stack(ims) - cell_roi = np.nanmax(ims, axis=0).astype(float) - cell_roi[cell_roi == 0] = np.nan - cell_roi -= 1 + if "dview" in locals(): + stop_server(dview=dview) # noqa: F821 - if cnm.estimates.b is not None and cnm.estimates.b.size != 0: - non_cell_roi_ims = get_roi( - scipy.sparse.csc_matrix(cnm.estimates.b), thr, thr_method, swap_dim, dims + c, dview, n_processes = setup_cluster( + backend="local", n_processes=None, single_thread=True ) - non_cell_roi_ims = np.stack(non_cell_roi_ims) - non_cell_roi = np.nanmax(non_cell_roi_ims, axis=0).astype(float) - else: - non_cell_roi_ims = None - non_cell_roi = np.zeros(dims) - non_cell_roi[non_cell_roi == 0] = np.nan - - all_roi = np.nanmax(np.stack([cell_roi, non_cell_roi]), axis=0) - - # NWBの追加 - nwbfile = {} - # NWBにROIを追加 - roi_list = [] - n_cells = cnm.estimates.A.shape[-1] - for i in range(n_cells): - kargs = {} - kargs["image_mask"] = cnm.estimates.A.T[i].T.toarray().reshape(dims) - if hasattr(cnm.estimates, "accepted_list"): - kargs["accepted"] = i in cnm.estimates.accepted_list - if hasattr(cnm.estimates, "rejected_list"): - kargs["rejected"] = i in cnm.estimates.rejected_list - roi_list.append(kargs) - - # backgroundsを追加 - if cnm.estimates.b is not None: - for bg in cnm.estimates.b.T: + + cnm = cnmf.CNMF(n_processes=n_processes, dview=dview, Ain=Ain, params=ops) + cnm = cnm.fit(mmap_images) + + if do_refit: + cnm = cnm.refit(mmap_images, dview=dview) + + stop_server(dview=dview) + + # contours plot + Cn = local_correlations(mmap_images.transpose(1, 2, 0)) + Cn[np.isnan(Cn)] = 0 + + thr_method = "nrg" + swap_dim = False + + iscell = np.concatenate( + [ + np.ones(cnm.estimates.A.shape[-1]), + np.zeros( + cnm.estimates.b.shape[-1] if cnm.estimates.b is not None else 0 + ), + ] + ).astype(bool) + + ims = self.get_roi(cnm.estimates.A, thr, thr_method, swap_dim, dims) + ims = np.stack(ims) + cell_roi = np.nanmax(ims, axis=0).astype(float) + cell_roi[cell_roi == 0] = np.nan + cell_roi -= 1 + + if cnm.estimates.b is not None and cnm.estimates.b.size != 0: + non_cell_roi_ims = self.get_roi( + scipy.sparse.csc_matrix(cnm.estimates.b), + thr, + thr_method, + swap_dim, + dims, + ) + non_cell_roi_ims = np.stack(non_cell_roi_ims) + non_cell_roi = np.nanmax(non_cell_roi_ims, axis=0).astype(float) + else: + non_cell_roi_ims = None + non_cell_roi = np.zeros(dims) + non_cell_roi[non_cell_roi == 0] = np.nan + + all_roi = np.nanmax(np.stack([cell_roi, non_cell_roi]), axis=0) + + # NWBの追加 + nwbfile = {} + # NWBにROIを追加 + roi_list = [] + n_cells = cnm.estimates.A.shape[-1] + for i in range(n_cells): kargs = {} - kargs["image_mask"] = bg.reshape(dims) + kargs["image_mask"] = cnm.estimates.A.T[i].T.toarray().reshape(dims) if hasattr(cnm.estimates, "accepted_list"): - kargs["accepted"] = False + kargs["accepted"] = i in cnm.estimates.accepted_list if hasattr(cnm.estimates, "rejected_list"): - kargs["rejected"] = False + kargs["rejected"] = i in cnm.estimates.rejected_list roi_list.append(kargs) - nwbfile[NWBDATASET.ROI] = {function_id: roi_list} - - # iscellを追加 - nwbfile[NWBDATASET.COLUMN] = { - function_id: { - "name": "iscell", - "description": "two columns - iscell & probcell", - "data": iscell, + # backgroundsを追加 + if cnm.estimates.b is not None: + for bg in cnm.estimates.b.T: + kargs = {} + kargs["image_mask"] = bg.reshape(dims) + if hasattr(cnm.estimates, "accepted_list"): + kargs["accepted"] = False + if hasattr(cnm.estimates, "rejected_list"): + kargs["rejected"] = False + roi_list.append(kargs) + + nwbfile[NWBDATASET.ROI] = {self.function_id: roi_list} + + # iscellを追加 + nwbfile[NWBDATASET.COLUMN] = { + self.function_id: { + "name": "iscell", + "description": "two columns - iscell & probcell", + "data": iscell, + } } - } - # Fluorescence - n_rois = len(cnm.estimates.C) - n_bg = len(cnm.estimates.f) if cnm.estimates.f is not None else 0 + # Fluorescence + n_rois = len(cnm.estimates.C) + n_bg = len(cnm.estimates.f) if cnm.estimates.f is not None else 0 - fluorescence = ( - np.concatenate( - [ - cnm.estimates.C, - cnm.estimates.f, - ] + fluorescence = ( + np.concatenate( + [ + cnm.estimates.C, + cnm.estimates.f, + ] + ) + if cnm.estimates.f is not None + else cnm.estimates.C ) - if cnm.estimates.f is not None - else cnm.estimates.C - ) - nwbfile[NWBDATASET.FLUORESCENCE] = { - function_id: { - "Fluorescence": { - "table_name": "ROIs", - "region": list(range(n_rois + n_bg)), - "name": "Fluorescence", - "data": fluorescence.T, - "unit": "lumens", + nwbfile[NWBDATASET.FLUORESCENCE] = { + self.function_id: { + "Fluorescence": { + "table_name": "ROIs", + "region": list(range(n_rois + n_bg)), + "name": "Fluorescence", + "data": fluorescence.T, + "unit": "lumens", + } } } - } - - cnmf_data = {} - cnmf_data["fluorescence"] = fluorescence - cnmf_data["im"] = ( - np.concatenate([ims, non_cell_roi_ims], axis=0) - if non_cell_roi_ims is not None - else ims - ) - cnmf_data["is_cell"] = iscell.astype(bool) - cnmf_data["images"] = mmap_images - - info = { - "images": ImageData( - np.array(Cn * 255, dtype=np.uint8), - output_dir=output_dir, - file_name="images", - ), - "fluorescence": FluoData(fluorescence, file_name="fluorescence"), - "iscell": IscellData(iscell, file_name="iscell"), - "all_roi": RoiData(all_roi, output_dir=output_dir, file_name="all_roi"), - "cell_roi": RoiData(cell_roi, output_dir=output_dir, file_name="cell_roi"), - "non_cell_roi": RoiData( - non_cell_roi, output_dir=output_dir, file_name="non_cell_roi" - ), - "nwbfile": nwbfile, - "cnmf_data": CaimanCnmfData(cnmf_data), - } - return info + cnmf_data = {} + cnmf_data["fluorescence"] = fluorescence + cnmf_data["im"] = ( + np.concatenate([ims, non_cell_roi_ims], axis=0) + if non_cell_roi_ims is not None + else ims + ) + cnmf_data["is_cell"] = iscell.astype(bool) + cnmf_data["images"] = mmap_images + + info = { + "images": ImageData( + np.array(Cn * 255, dtype=np.uint8), + output_dir=self.output_dir, + file_name="images", + ), + "fluorescence": FluoData(fluorescence, file_name="fluorescence"), + "iscell": IscellData(iscell, file_name="iscell"), + "all_roi": RoiData( + all_roi, output_dir=self.output_dir, file_name="all_roi" + ), + "cell_roi": RoiData( + cell_roi, output_dir=self.output_dir, file_name="cell_roi" + ), + "non_cell_roi": RoiData( + non_cell_roi, output_dir=self.output_dir, file_name="non_cell_roi" + ), + "nwbfile": nwbfile, + "cnmf_data": CaimanCnmfData(cnmf_data), + } + + return info + + @staticmethod + def get_roi(A, thr, thr_method, swap_dim, dims): + from scipy.ndimage import binary_fill_holes + from skimage.measure import find_contours + + d, nr = np.shape(A) + + # for each patches + ims = [] + coordinates = [] + for i in range(nr): + pars = dict() + # we compute the cumulative sum of the energy of the Ath component + # that has been ordered from least to highest + patch_data = A.data[A.indptr[i] : A.indptr[i + 1]] + indx = np.argsort(patch_data)[::-1] + + if thr_method == "nrg": + cumEn = np.cumsum(patch_data[indx] ** 2) + if len(cumEn) == 0: + pars = dict( + coordinates=np.array([]), + CoM=np.array([np.NaN, np.NaN]), + neuron_id=i + 1, + ) + coordinates.append(pars) + continue + else: + # we work with normalized values + cumEn /= cumEn[-1] + Bvec = np.ones(d) + # we put it in a similar matrix + Bvec[A.indices[A.indptr[i] : A.indptr[i + 1]][indx]] = cumEn + else: + Bvec = np.zeros(d) + Bvec[A.indices[A.indptr[i] : A.indptr[i + 1]]] = ( + patch_data / patch_data.max() + ) + + if swap_dim: + Bmat = np.reshape(Bvec, dims, order="C") + else: + Bmat = np.reshape(Bvec, dims, order="F") + + r_mask = np.zeros_like(Bmat, dtype="bool") + contour = find_contours(Bmat, thr) + for c in contour: + r_mask[ + np.round(c[:, 0]).astype("int"), np.round(c[:, 1]).astype("int") + ] = 1 + + # Fill in the hole created by the contour boundary + r_mask = binary_fill_holes(r_mask) + ims.append(r_mask + (i * r_mask)) + + return ims diff --git a/studio/app/optinist/wrappers/caiman/cnmfe.py b/studio/app/optinist/wrappers/caiman/cnmfe.py index 60cf4dac4..0f47b1854 100644 --- a/studio/app/optinist/wrappers/caiman/cnmfe.py +++ b/studio/app/optinist/wrappers/caiman/cnmfe.py @@ -1,17 +1,141 @@ +from typing import List, Optional + +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import ImageData from studio.app.optinist.dataclass import FluoData, IscellData -from studio.app.optinist.wrappers.caiman.cnmf import caiman_cnmf +from studio.app.optinist.wrappers.caiman.cnmf import CaimanCnmf + + +class CaimanCnmfEInitParams(BaseModel): + # Ain: Optional[List] = Field( + # None, description="possibility to seed with predetermined binary masks" + # ) + do_refit: bool = Field( + False, + description=( + "Whether to re-run seeded CNMFE on accepted patches " + "to refine and perform deconvolution" + ), + ) + K: Optional[int] = Field( + 4, description="upper bound on number of components per patch, in general None" + ) + gSig: List[int] = Field( + [3, 3], + description=( + "gaussian width of a 2D gaussian kernel, which approximates a neuron" + ), + ) + gSiz: List[int] = Field( + [13, 13], + description=("average diameter of a neuron, in general 4*gSig+1"), + ) + ssub: int = Field( + 1, + description=( + "downsampling factor in space for initialization, " + "increase if you have memory problems. " + "you can pass them here as boolean vectors" + ), + ) + tsub: int = Field( + 1, + description=( + "downsampling factor in time for initialization, " + "increase if you have memory problems." + ), + ) + nb: int = Field( + 2, + description=""" + number of background components (rank) if positive, + else exact ring model with following settings. + + gnb= 0: Return background as b and W + + gnb=-1: Return full rank background B + + gnb<-1: Don't return background + """, + ) + min_corr: float = Field(0.8, description="min peak value from correlation image") + min_pnr: int = Field(10, description="min peak to noise ration from PNR image") + ssub_B: int = Field( + 2, description="additional downsampling factor in space for background" + ) + ring_size_factor: float = Field( + 1.4, description=("radius of ring is gSiz*ring_size_factor") + ) + + +class CaimanCnmfEPreprocessParams(BaseModel): + p: int = Field(1, description="order of the autoregressive system") + + +class CaimanCnmfEPatchParams(BaseModel): + rf: Optional[int] = Field( + 40, + description=( + "half-size of the patches in pixels. e.g., if rf=40, patches are 80x80" + ), + ) + stride: int = Field( + 20, + description=( + "amount of overlap between the patches in pixels, " + "(keep it at least large as gSiz, i.e 4 times the neuron size gSig)" + ), + ) + low_rank_background: Optional[bool] = Field( + None, + description=( + "None leaves background of each patch intact, " + "True performs global low-rank approximation if gnb>0" + ), + ) + nb_patch: int = Field( + 0, + description=( + "number of background components (rank) per patch if gnb>0, " + "else it is set automatically" + ), + ) + + +class CaimanCnmfEMergeParams(BaseModel): + thr: float = Field(0.9) + merge_thr: float = Field( + 0.7, description="merging threshold, max correlation allowed" + ) + + +class CaimanCnmfEParams(BaseModel): + init_params: CaimanCnmfEInitParams = Field(CaimanCnmfEInitParams()) + preprocess_params: CaimanCnmfEPreprocessParams = Field( + CaimanCnmfEPreprocessParams() + ) + patch_params: CaimanCnmfEPatchParams = Field(CaimanCnmfEPatchParams()) + merge_params: CaimanCnmfEMergeParams = Field(CaimanCnmfEMergeParams()) -def caiman_cnmfe( - images: ImageData, output_dir: str, params: dict = None, **kwargs -) -> dict(fluorescence=FluoData, iscell=IscellData): - cnmfe_fixed_params = { - "center_psf": True, - "method_init": "corr_pnr", # use this for 1 photon - "only_init": True, # set it to True to run CNMF-E - "normalize_init": False, - } - params["fixed_params"] = cnmfe_fixed_params +class CaimanCnmfE(AlgoTemplate): + def run( + self, params: CaimanCnmfEParams, images: ImageData + ) -> dict(fluorescence=FluoData, iscell=IscellData): + cnmfe_fixed_params = { + "center_psf": True, + "method_init": "corr_pnr", # use this for 1 photon + "only_init": True, # set it to True to run CNMF-E + "normalize_init": False, + } + params["fixed_params"] = cnmfe_fixed_params - return caiman_cnmf(images=images, output_dir=output_dir, params=params, **kwargs) + return ( + CaimanCnmf() + .set_output_dir(self.output_dir) + .set_nwb_params(self.nwb_params) + .run(params=params, images=images) + ) diff --git a/studio/app/optinist/wrappers/caiman/motion_correction.py b/studio/app/optinist/wrappers/caiman/motion_correction.py index aa9c98692..ee91f2c38 100644 --- a/studio/app/optinist/wrappers/caiman/motion_correction.py +++ b/studio/app/optinist/wrappers/caiman/motion_correction.py @@ -1,5 +1,10 @@ import shutil +from typing import List, Optional +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.core.utils.filepath_creater import ( create_directory, join_filepath, @@ -9,91 +14,115 @@ from studio.app.optinist.dataclass import RoiData -def caiman_mc( - image: ImageData, output_dir: str, params: dict = None, **kwargs -) -> dict(mc_images=ImageData): - import numpy as np - from caiman import load_memmap, save_memmap, stop_server - from caiman.base.rois import extract_binary_masks_from_structural_channel - from caiman.cluster import setup_cluster - from caiman.motion_correction import MotionCorrect - from caiman.source_extraction.cnmf.params import CNMFParams - - function_id = output_dir.split("/")[-1] - print("start caiman motion_correction:", function_id) - - opts = CNMFParams() - - if params is not None: - opts.change_params(params_dict=params) - - c, dview, n_processes = setup_cluster( - backend="local", n_processes=None, single_thread=True - ) - - mc = MotionCorrect(image.path, dview=dview, **opts.get_group("motion")) - - mc.motion_correct(save_movie=True) - border_to_0 = 0 if mc.border_nan == "copy" else mc.border_to_0 - - # memory mapping - fname_new = save_memmap( - mc.mmap_file, base_name=function_id, order="C", border_to_0=border_to_0 - ) - - stop_server(dview=dview) - - # now load the file - Yr, dims, T = load_memmap(fname_new) - - images = np.array(Yr.T.reshape((T,) + dims, order="F")) - - meanImg = images.mean(axis=0) - rois = ( - extract_binary_masks_from_structural_channel( - meanImg, gSig=7, expand_method="dilation" - )[0] - .reshape(meanImg.shape[0], meanImg.shape[1], -1) - .transpose(2, 0, 1) - ) - - rois = rois.astype(np.float) - - for i, _ in enumerate(rois): - rois[i] *= i + 1 - - rois = np.nanmax(rois, axis=0) - rois[rois == 0] = np.nan - rois -= 1 +class CaimanMcParams(BaseModel): + border_nan: str = Field("copy") + gSig_filt: Optional[List] = Field(None) + is3D: bool = Field(False) + max_deviation_rigid: int = Field(3) + max_shifts: List[int] = Field([6, 6]) + min_mov: Optional[float] = Field(None) + niter_rig: int = Field(1) + nonneg_movie: bool = Field(True) + num_frames_split: int = Field(80) + num_splits_to_process_els: Optional[int] = Field(None) + num_splits_to_process_rig: Optional[int] = Field(None) + overlaps: List[int] = Field([32, 32]) + pw_rigid: bool = Field(False) + shifts_opencv: bool = Field(True) + splits_els: int = Field(14) + splits_rig: int = Field(14) + strides: List[int] = Field([96, 96]) + upsample_factor_grid: int = Field(4) + use_cuda: bool = Field(False) + + +class CaimanMc(AlgoTemplate): + def run( + self, params: CaimanMcParams, image: ImageData + ) -> dict(mc_images=ImageData): + import numpy as np + from caiman import load_memmap, save_memmap, stop_server + from caiman.base.rois import extract_binary_masks_from_structural_channel + from caiman.cluster import setup_cluster + from caiman.motion_correction import MotionCorrect + from caiman.source_extraction.cnmf.params import CNMFParams + + print("start caiman motion_correction:", self.function_id) + + opts = CNMFParams() + + if params is not None: + opts.change_params(params_dict=params) + + c, dview, n_processes = setup_cluster( + backend="local", n_processes=None, single_thread=True + ) + + mc = MotionCorrect(image.path, dview=dview, **opts.get_group("motion")) + + mc.motion_correct(save_movie=True) + border_to_0 = 0 if mc.border_nan == "copy" else mc.border_to_0 + + # memory mapping + fname_new = save_memmap( + mc.mmap_file, base_name=self.function_id, order="C", border_to_0=border_to_0 + ) + + stop_server(dview=dview) + + # now load the file + Yr, dims, T = load_memmap(fname_new) + + images = np.array(Yr.T.reshape((T,) + dims, order="F")) + + meanImg = images.mean(axis=0) + rois = ( + extract_binary_masks_from_structural_channel( + meanImg, gSig=7, expand_method="dilation" + )[0] + .reshape(meanImg.shape[0], meanImg.shape[1], -1) + .transpose(2, 0, 1) + ) + + rois = rois.astype(np.float) + + for i, _ in enumerate(rois): + rois[i] *= i + 1 + + rois = np.nanmax(rois, axis=0) + rois[rois == 0] = np.nan + rois -= 1 + + xy_trans_data = ( + (np.array(mc.x_shifts_els), np.array(mc.y_shifts_els)) + if params["pw_rigid"] + else np.array(mc.shifts_rig) + ) + + mc_images = ImageData(images, output_dir=self.output_dir, file_name="mc_images") + + nwbfile = {} + nwbfile[NWBDATASET.MOTION_CORRECTION] = { + self.function_id: { + "mc_data": mc_images, + "xy_trans_data": xy_trans_data, + } + } - xy_trans_data = ( - (np.array(mc.x_shifts_els), np.array(mc.y_shifts_els)) - if params["pw_rigid"] - else np.array(mc.shifts_rig) - ) + info = { + "mc_images": mc_images, + "meanImg": ImageData( + meanImg, output_dir=self.output_dir, file_name="meanImg" + ), + "rois": RoiData(rois, output_dir=self.output_dir, file_name="rois"), + "nwbfile": nwbfile, + } - mc_images = ImageData(images, output_dir=output_dir, file_name="mc_images") + # Clean up temporary files + mmap_output_dir = join_filepath([self.output_dir, "mmap"]) + create_directory(mmap_output_dir) + for mmap_file in mc.mmap_file: + shutil.move(mmap_file, mmap_output_dir) + shutil.move(fname_new, mmap_output_dir) - nwbfile = {} - nwbfile[NWBDATASET.MOTION_CORRECTION] = { - function_id: { - "mc_data": mc_images, - "xy_trans_data": xy_trans_data, - } - } - - info = { - "mc_images": mc_images, - "meanImg": ImageData(meanImg, output_dir=output_dir, file_name="meanImg"), - "rois": RoiData(rois, output_dir=output_dir, file_name="rois"), - "nwbfile": nwbfile, - } - - # Clean up temporary files - mmap_output_dir = join_filepath([output_dir, "mmap"]) - create_directory(mmap_output_dir) - for mmap_file in mc.mmap_file: - shutil.move(mmap_file, mmap_output_dir) - shutil.move(fname_new, mmap_output_dir) - - return info + return info diff --git a/studio/app/optinist/wrappers/caiman/params/caiman_cnmf.yaml b/studio/app/optinist/wrappers/caiman/params/caiman_cnmf.yaml deleted file mode 100644 index 4803d65e7..000000000 --- a/studio/app/optinist/wrappers/caiman/params/caiman_cnmf.yaml +++ /dev/null @@ -1,21 +0,0 @@ -init_params: - # Ain: null # TBD: need to support 2D array type. - do_refit: False - - K: 4 - gSig: [4, 4] - ssub: 1 - tsub: 1 - nb: 2 - method_init: "greedy_roi" - -preprocess_params: - p: 1 - -patch_params: - rf: - stride: 6 - -merge_params: - thr: 0.9 - merge_thr: 0.85 diff --git a/studio/app/optinist/wrappers/caiman/params/caiman_cnmfe.yaml b/studio/app/optinist/wrappers/caiman/params/caiman_cnmfe.yaml deleted file mode 100644 index 199f25923..000000000 --- a/studio/app/optinist/wrappers/caiman/params/caiman_cnmfe.yaml +++ /dev/null @@ -1,38 +0,0 @@ -init_params: - # Ain: null # possibility to seed with predetermined binary masks - # TBD: need to support 2D array type. - do_refit: False - - K: 4 # upper bound on number of components per patch, in general None - gSig: [3, 3] # gaussian width of a 2D gaussian kernel, which approximates a neuron - gSiz: [13, 13] # average diameter of a neuron, in general 4*gSig+1 - ssub: 1 # downsampling factor in space for initialization, - # increase if you have memory problems - # you can pass them here as boolean vectors - tsub: 2 # downsampling factor in time for initialization, - # increase if you have memory problems - nb: 0 # number of background components (rank) if positive, - # else exact ring model with following settings - # gnb= 0: Return background as b and W - # gnb=-1: Return full rank background B - # gnb<-1: Don't return background - min_corr: .8 # min peak value from correlation image - min_pnr: 10 # min peak to noise ration from PNR image - ssub_B: 2 # additional downsampling factor in space for background - ring_size_factor: 1.4 # radius of ring is gSiz*ring_size_factor - -preprocess_params: - p: 1 # order of the autoregressive system - -patch_params: - rf: 40 # half-size of the patches in pixels. e.g., if rf=40, patches are 80x80 - stride: 20 # amount of overlap between the patches in pixels - # (keep it at least large as gSiz, i.e 4 times the neuron size gSig) - low_rank_background: null # None leaves background of each patch intact, - # True performs global low-rank approximation if gnb>0 - nb_patch: 0 # number of background components (rank) per patch if gnb>0, - # else it is set automatically - -merge_params: - thr: 0.9 - merge_thr: .7 # merging threshold, max correlation allowed diff --git a/studio/app/optinist/wrappers/caiman/params/caiman_mc.yaml b/studio/app/optinist/wrappers/caiman/params/caiman_mc.yaml deleted file mode 100644 index 6f6eb462f..000000000 --- a/studio/app/optinist/wrappers/caiman/params/caiman_mc.yaml +++ /dev/null @@ -1,19 +0,0 @@ -border_nan: 'copy' -gSig_filt: null -is3D: False -max_deviation_rigid: 3 -max_shifts: [6, 6] -min_mov: null -niter_rig: 1 -nonneg_movie: True -num_frames_split: 80 -num_splits_to_process_els: null -num_splits_to_process_rig: null -overlaps: [32, 32] -pw_rigid: False -shifts_opencv: True -splits_els: 14 -splits_rig: 14 -strides: [96, 96] -upsample_factor_grid: 4 -use_cuda: False diff --git a/studio/app/optinist/wrappers/lccd/__init__.py b/studio/app/optinist/wrappers/lccd/__init__.py index c9152dd0c..1f644b55d 100644 --- a/studio/app/optinist/wrappers/lccd/__init__.py +++ b/studio/app/optinist/wrappers/lccd/__init__.py @@ -1,9 +1,9 @@ -from studio.app.optinist.wrappers.lccd.lccd_detection import lccd_detect +from studio.app.optinist.wrappers.lccd.lccd_detection import LccdDetect lccd_wrapper_dict = { "lccd": { "lccd_cell_detection": { - "function": lccd_detect, + "function": LccdDetect, "conda_name": "lccd", }, } diff --git a/studio/app/optinist/wrappers/lccd/lccd_detection.py b/studio/app/optinist/wrappers/lccd/lccd_detection.py index ddaa712d6..461178475 100644 --- a/studio/app/optinist/wrappers/lccd/lccd_detection.py +++ b/studio/app/optinist/wrappers/lccd/lccd_detection.py @@ -1,108 +1,148 @@ import numpy as np +from pydantic import BaseModel +from pydantic.dataclasses import Field +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import ImageData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import FluoData, LccdData, RoiData -def lccd_detect( - mc_images: ImageData, output_dir: str, params: dict = None, **kwargs -) -> dict(fluorescence=FluoData, cell_roi=RoiData): - from studio.app.optinist.wrappers.lccd.lccd_python.lccd import LCCD - - function_id = output_dir.split("/")[-1] - print("start lccd_detect:", function_id) - - print("params: ", params) - lccd = LCCD(params) - D = LoadData(mc_images) - assert len(D.shape) == 3, "input array should have dimensions (width, height, time)" - roi = lccd.apply(D) - - dff_f0_frames = params["dff"]["f0_frames"] - dff_f0_percentile = params["dff"]["f0_percentile"] - num_cell = roi.shape[1] - num_frames = D.shape[2] - is_cell = np.ones(num_cell, dtype=bool) - - reshapedD = D.reshape([D.shape[0] * D.shape[1], D.shape[2]]) - timeseries = np.zeros([num_cell, num_frames]) - roi_list = [] - - for i in range(num_cell): - roi_list.append((roi[:, i].reshape([D.shape[0], D.shape[1]])) * (i + 1)) - timeseries[i, :] = np.mean(reshapedD[roi[:, i] > 0, :], axis=0) - - im = np.stack(roi_list) - im[im == 0] = np.nan - im -= 1 - - timeseries_dff = np.ones([num_cell, num_frames]) * np.nan - for i in range(num_cell): - for k in range(num_frames): - if (k - dff_f0_frames >= 0) and (k + dff_f0_frames < num_frames): - f0 = np.percentile( - timeseries[i, k - dff_f0_frames : k + dff_f0_frames], - dff_f0_percentile, - ) - timeseries_dff[i, k] = (timeseries[i, k] - f0) / f0 - - roi_list = [{"image_mask": roi[:, i].reshape(D.shape[:2])} for i in range(num_cell)] - - nwbfile = {} - nwbfile[NWBDATASET.ROI] = {function_id: roi_list} - nwbfile[NWBDATASET.COLUMN] = { - function_id: { - "name": "iscell", - "description": "two columns - iscell & probcell", - "data": is_cell, - } - } - - nwbfile[NWBDATASET.FLUORESCENCE] = { - function_id: { - "Fluorescence": { - "table_name": "Fluorescence", - "region": list(range(len(timeseries))), - "name": "Fluorescence", - "data": timeseries, - "unit": "lumens", +class LccdBlobDetectorParams(BaseModel): + filtersize1: int = Field(100) + filtersize2: int = Field(4) + sigma: float = Field(1.25) + fsize: int = Field(30) + min_area: int = Field(20) + max_area: int = Field(50) + sparse: bool = Field(False) + + +class LccdRoiIntegrationParams(BaseModel): + overlap_threshold: float = Field(0.4) + min_area: int = Field(20) + max_area: int = Field(100) + sparse: bool = Field(False) + + +class LccdMainParams(BaseModel): + frame_divider: int = Field(100) + + +class LccdDffParams(BaseModel): + f0_frames: int = Field(100) + f0_percentile: int = Field(8) + + +class LccdDetectParams(BaseModel): + blob_detector: LccdBlobDetectorParams = Field(LccdBlobDetectorParams()) + roi_integration: LccdRoiIntegrationParams = Field(LccdRoiIntegrationParams()) + lccd: LccdMainParams = Field(LccdMainParams()) + dff: LccdDffParams = Field(LccdDffParams()) + + +class LccdDetect(AlgoTemplate): + def run( + self, params: LccdDetectParams, mc_images: ImageData + ) -> dict(fluorescence=FluoData, cell_roi=RoiData): + from studio.app.optinist.wrappers.lccd.lccd_python.lccd import LCCD + + print("start lccd_detect:", self.function_id) + + print("params: ", params) + lccd = LCCD(params) + D = self.LoadData(mc_images) + assert ( + len(D.shape) == 3 + ), "input array should have dimensions (width, height, time)" + roi = lccd.apply(D) + + dff_f0_frames = params["dff"]["f0_frames"] + dff_f0_percentile = params["dff"]["f0_percentile"] + num_cell = roi.shape[1] + num_frames = D.shape[2] + is_cell = np.ones(num_cell, dtype=bool) + + reshapedD = D.reshape([D.shape[0] * D.shape[1], D.shape[2]]) + timeseries = np.zeros([num_cell, num_frames]) + roi_list = [] + + for i in range(num_cell): + roi_list.append((roi[:, i].reshape([D.shape[0], D.shape[1]])) * (i + 1)) + timeseries[i, :] = np.mean(reshapedD[roi[:, i] > 0, :], axis=0) + + im = np.stack(roi_list) + im[im == 0] = np.nan + im -= 1 + + timeseries_dff = np.ones([num_cell, num_frames]) * np.nan + for i in range(num_cell): + for k in range(num_frames): + if (k - dff_f0_frames >= 0) and (k + dff_f0_frames < num_frames): + f0 = np.percentile( + timeseries[i, k - dff_f0_frames : k + dff_f0_frames], + dff_f0_percentile, + ) + timeseries_dff[i, k] = (timeseries[i, k] - f0) / f0 + + roi_list = [ + {"image_mask": roi[:, i].reshape(D.shape[:2])} for i in range(num_cell) + ] + + nwbfile = {} + nwbfile[NWBDATASET.ROI] = {self.function_id: roi_list} + nwbfile[NWBDATASET.COLUMN] = { + self.function_id: { + "name": "iscell", + "description": "two columns - iscell & probcell", + "data": is_cell, } } - } - - lccd_data = {} - lccd_data["images"] = D - lccd_data["roi"] = roi - lccd_data["is_cell"] = is_cell - info = { - "lccd": LccdData(lccd_data), - "cell_roi": RoiData( - np.nanmax(im, axis=0), output_dir=output_dir, file_name="cell_roi" - ), - "fluorescence": FluoData(timeseries, file_name="fluorescence"), - "dff": FluoData(timeseries_dff, file_name="dff"), - "nwbfile": nwbfile, - } + nwbfile[NWBDATASET.FLUORESCENCE] = { + self.function_id: { + "Fluorescence": { + "table_name": "Fluorescence", + "region": list(range(len(timeseries))), + "name": "Fluorescence", + "data": timeseries, + "unit": "lumens", + } + } + } - return info + lccd_data = {} + lccd_data["images"] = D + lccd_data["roi"] = roi + lccd_data["is_cell"] = is_cell + + info = { + "lccd": LccdData(lccd_data), + "cell_roi": RoiData( + np.nanmax(im, axis=0), output_dir=self.output_dir, file_name="cell_roi" + ), + "fluorescence": FluoData(timeseries, file_name="fluorescence"), + "dff": FluoData(timeseries_dff, file_name="dff"), + "nwbfile": nwbfile, + } + return info -def LoadData(mc_images): - from PIL import Image + @staticmethod + def LoadData(mc_images): + from PIL import Image - img_pile = Image.open(mc_images.path[0]) - num_page = img_pile.n_frames + img_pile = Image.open(mc_images.path[0]) + num_page = img_pile.n_frames - imgs = np.zeros((img_pile.size[0], img_pile.size[1], num_page)) - for i in range(num_page): - img_pile.seek(i) - imgs[:, :, i] = np.asarray(img_pile) - img_pile.close() + imgs = np.zeros((img_pile.size[0], img_pile.size[1], num_page)) + for i in range(num_page): + img_pile.seek(i) + imgs[:, :, i] = np.asarray(img_pile) + img_pile.close() - maxval = np.max(imgs) - minval = np.min(imgs) - imgs = (imgs - minval) / (maxval - minval) + maxval = np.max(imgs) + minval = np.min(imgs) + imgs = (imgs - minval) / (maxval - minval) - return imgs + return imgs diff --git a/studio/app/optinist/wrappers/lccd/lccd_python/roi_integration.py b/studio/app/optinist/wrappers/lccd/lccd_python/roi_integration.py index 60b869701..076eb2a02 100644 --- a/studio/app/optinist/wrappers/lccd/lccd_python/roi_integration.py +++ b/studio/app/optinist/wrappers/lccd/lccd_python/roi_integration.py @@ -6,18 +6,17 @@ def filter_roi_by_area(roi, min_area, max_area): """ - Parameters - ---------- - roi: np.array or scipy.sparse.csc (Maybe work with other sparse matrix types) - roi has shape (X, n). - X is spatial dimension. - n is the number of rois - min_area: int - max_area: int - - Returns - ------- - roi + Args: + roi (np.array or scipy.sparse.csc): Maybe work with other sparse matrix types. + roi has shape (X, n). + X is spatial dimension. + n is the number of rois + min_area (int) + max_area (int) + + Returns: + roi + """ n_rois = roi.shape[1] @@ -80,14 +79,17 @@ def remove_overlap(self, roi_a, roi_b, i, j): From both roi_a's ith region and roi_b's jth region, remove their intersection. - Let - roi_a[:, i] = [1, 0, 1, 1, 0, 1, 1].T - roi_b[:, j] = [0, 1, 1, 0, 1, 1, 0].T + :: + + Let + roi_a[:, i] = [1, 0, 1, 1, 0, 1, 1].T + roi_b[:, j] = [0, 1, 1, 0, 1, 1, 0].T + + intersection is [0, 0, 1, 0, 0, 1, 0].T + then, calling remove_overlap(roi_a, roi_b, i, j), + roi_a[:, i] = [1, 0, 0, 1, 0, 0, 1].T + roi_b[:, j] = [0, 1, 0, 0, 1, 0, 0].T - intersection is [0, 0, 1, 0, 0, 1, 0].T - then, calling remove_overlap(roi_a, roi_b, i, j), - roi_a[:, i] = [1, 0, 0, 1, 0, 0, 1].T - roi_b[:, j] = [0, 1, 0, 0, 1, 0, 0].T """ if self.sparse: # TODO: improve readability diff --git a/studio/app/optinist/wrappers/lccd/lccd_python/utils.py b/studio/app/optinist/wrappers/lccd/lccd_python/utils.py index 3248823d7..f1123432e 100644 --- a/studio/app/optinist/wrappers/lccd/lccd_python/utils.py +++ b/studio/app/optinist/wrappers/lccd/lccd_python/utils.py @@ -25,47 +25,48 @@ def insert_binary_col_binary_csc(mat, indices, i=None): This method is in-place operation. - Parameters - ---------- - mat: scipy.sparse.csc_matrix - i: int - indices: numpy.array - The indices of nonzero elements in the column to be inserted. - - - Let mat = scipy.sparse.csc_matrix(np.array([ - [1, 0, 1], - [1, 1, 1], - [0, 0, 0], - ])) - i = 1 - indices = np.array([2]). - The column to be inserted is - [ - [0], - [0], - [1], - ] - insert_binary_col_binary_csc(mat, indices, i) changes mat to be - [ - [1, 0, 0, 1], - [1, 0, 1, 1], - [0, 1, 0, 0], - ] ^ Inserted column - - When indices = np.array([0, 2]), - The column to be inserted is - [ - [1], - [0], - [1], - ] - and insert_binary_col_binary_csc(mat, indices, i) changes mat to be - [ - [1, 1, 0, 1], - [1, 0, 1, 1], - [0, 1, 0, 0], - ] ^ Inserted column + :: + + Let mat = scipy.sparse.csc_matrix(np.array([ + [1, 0, 1], + [1, 1, 1], + [0, 0, 0], + ])) + i = 1 + indices = np.array([2]). + The column to be inserted is + [ + [0], + [0], + [1], + ] + insert_binary_col_binary_csc(mat, indices, i) changes mat to be + [ + [1, 0, 0, 1], + [1, 0, 1, 1], + [0, 1, 0, 0], + ] ^ Inserted column + + When indices = np.array([0, 2]), + The column to be inserted is + [ + [1], + [0], + [1], + ] + and insert_binary_col_binary_csc(mat, indices, i) changes mat to be + [ + [1, 1, 0, 1], + [1, 0, 1, 1], + [0, 1, 0, 0], + ] ^ Inserted column + + Args: + mat (scipy.sparse.csc_matrix) + i (int) + indices (numpy.array): The indices of nonzero elements in the column to be + inserted. + """ if not isinstance(mat, scipy.sparse.csc_matrix): @@ -191,50 +192,49 @@ def array_to_lil_row(arr): lil style sparse matrix of the large matrix is achieved. - Parameters - ---------- - arr: 1 dim numpy.array - - Returns - ------- - datum: list of lists - row: list of lists - - Examples - -------- - >>> arr = np.array([1, 2, 3]) - >>> array_row_to_lil_row(arr) - list([1, 2, 3]), list([0, 1, 2]) - - >>> arr = np.array([0, 0, 1]) - >>> array_row_to_lil_row(arr) - list([1]), list([2]) - - >>> arr = np.array([0, 0, 0]) - >>> array_row_to_lil_row(arr) - list([]), list([]) - - # example of array to sparse.lil_matrix conversion. - >>> arr = np.array([ - ... [1, 2, 3], - ... [0, 0, 1], - ... [0, 0, 0], - ... [7, 8, 0], - ... ]) - >>> rows, data = [], [] - >>> for row_ in arr: - ... datum, row = array_to_lil_row(row_) - ... rows.append(row) - ... data.append(datum) - ... n_rows, n_cols = 4, 3 - ... lil_arr = scipy.sparse.lil_matrix((n_rows, n_cols), dtype=data[0][0].dtype) - ... lil_arr.data = np.array(data, dtype='object') - ... lil_arr.rows = np.array(rows, dtype='object') - >>> lil_arr.toarray() - array([[1, 2, 3], - [0, 0, 1], - [0, 0, 0], - [7, 8, 0]]) + Args: + arr (numpy.array): 1-dim array + + Returns: + datum: list of lists + row: list of lists + + Examples:: + + >>> arr = np.array([1, 2, 3]) + >>> array_row_to_lil_row(arr) + list([1, 2, 3]), list([0, 1, 2]) + + >>> arr = np.array([0, 0, 1]) + >>> array_row_to_lil_row(arr) + list([1]), list([2]) + + >>> arr = np.array([0, 0, 0]) + >>> array_row_to_lil_row(arr) + list([]), list([]) + + # example of array to sparse.lil_matrix conversion. + >>> arr = np.array([ + ... [1, 2, 3], + ... [0, 0, 1], + ... [0, 0, 0], + ... [7, 8, 0], + ... ]) + >>> rows, data = [], [] + >>> for row_ in arr: + ... datum, row = array_to_lil_row(row_) + ... rows.append(row) + ... data.append(datum) + ... n_rows, n_cols = 4, 3 + ... lil_arr = scipy.sparse.lil_matrix((n_rows, n_cols), dtype=data[0][0].dtype) + ... lil_arr.data = np.array(data, dtype='object') + ... lil_arr.rows = np.array(rows, dtype='object') + >>> lil_arr.toarray() + array([[1, 2, 3], + [0, 0, 1], + [0, 0, 0], + [7, 8, 0]]) + """ row = list(np.where(arr != 0)[0]) diff --git a/studio/app/optinist/wrappers/lccd/params/lccd_cell_detection.yaml b/studio/app/optinist/wrappers/lccd/params/lccd_cell_detection.yaml deleted file mode 100644 index e0de10c48..000000000 --- a/studio/app/optinist/wrappers/lccd/params/lccd_cell_detection.yaml +++ /dev/null @@ -1,18 +0,0 @@ -blob_detector: - filtersize1: 100 - filtersize2: 4 - sigma: 1.25 - fsize: 30 - min_area: 20 - max_area: 50 - sparse: False -roi_integration: - overlap_threshold: 0.4 - min_area: 20 - max_area: 100 - sparse: False -lccd: - frame_divider: 100 -dff: - f0_frames: 100 - f0_percentile: 8 diff --git a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/__init__.py b/studio/app/optinist/wrappers/optinist/basic_neural_analysis/__init__.py index fc7398e89..102cf5063 100644 --- a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/__init__.py +++ b/studio/app/optinist/wrappers/optinist/basic_neural_analysis/__init__.py @@ -1,11 +1,9 @@ # from studio.app.optinist.wrappers.optinist.basic_neural_analysis.cell_grouping import ( # noqa: E501 -# cell_grouping, +# CellGrouping, # ) from studio.app.optinist.wrappers.optinist.basic_neural_analysis.eta import ETA basic_neural_analysis_wrapper_dict = { "eta": {"function": ETA}, - # 'cell_grouping': { - # 'function': cell_grouping - # } + # "cell_grouping": {"function": CellGrouping}, } diff --git a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/cell_grouping.py b/studio/app/optinist/wrappers/optinist/basic_neural_analysis/cell_grouping.py index e0622f842..bc2fb6eb1 100644 --- a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/cell_grouping.py +++ b/studio/app/optinist/wrappers/optinist/basic_neural_analysis/cell_grouping.py @@ -1,30 +1,35 @@ +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import TimeSeriesData -from studio.app.optinist.dataclass import NWBFile -def cell_grouping( - neural_data: TimeSeriesData, - output_dir: str, - nwbfile: NWBFile = None, - params: dict = None, - **kwargs, -) -> dict(): - import numpy as np +class CellGroupingParams(BaseModel): + transpose: bool = Field(False) + threshold: float = Field(1.0) + start_time: int = Field(-10) + end_time: int = Field(0) + + +class CellGrouping(AlgoTemplate): + def run(self, params: CellGroupingParams, neural_data: TimeSeriesData) -> dict(): + import numpy as np - neural_data = neural_data.data - std = neural_data.std + neural_data = neural_data.data + std = neural_data.std - if params["transpose"]: - neural_data = neural_data.transpose() + if params["transpose"]: + neural_data = neural_data.transpose() - baseline = np.mean(neural_data, axis=1, keepdims=True) - std = neural_data[np.argmax(neural_data, axis=1)].std() + baseline = np.mean(neural_data, axis=1, keepdims=True) + std = neural_data[np.argmax(neural_data, axis=1)].std() - grouped_cells = (neural_data.max(axis=1) - baseline) / std > params["threshold"] + grouped_cells = (neural_data.max(axis=1) - baseline) / std > params["threshold"] - info = {} - info["grouped_cells"] = TimeSeriesData( - neural_data[grouped_cells], std=std, file_name="grouped_cells_mean" - ) + info = {} + info["grouped_cells"] = TimeSeriesData( + neural_data[grouped_cells], std=std, file_name="grouped_cells_mean" + ) - return info + return info diff --git a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/eta.py b/studio/app/optinist/wrappers/optinist/basic_neural_analysis/eta.py index e3fbd318a..b57e14add 100644 --- a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/eta.py +++ b/studio/app/optinist/wrappers/optinist/basic_neural_analysis/eta.py @@ -1,124 +1,138 @@ import numpy as np +from pydantic import BaseModel +from pydantic.dataclasses import Field +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import HeatMapData, TimeSeriesData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import BehaviorData, FluoData, IscellData -def calc_trigger(behavior_data, trigger_type, trigger_threshold): - flg = np.array(behavior_data > trigger_threshold, dtype=int) - if trigger_type == "up": - trigger_idx = np.where(np.ediff1d(flg) == 1)[0] - elif trigger_type == "down": - trigger_idx = np.where(np.ediff1d(flg) == -1)[0] - elif trigger_type == "cross": - trigger_idx = np.where(np.ediff1d(flg) != 0)[0] - else: - trigger_idx = np.where(np.ediff1d(flg) == 0)[0] - - return trigger_idx - - -def calc_trigger_average(neural_data, trigger_idx, start_time, end_time): - num_frame = neural_data.shape[0] - - ind = np.array(range(start_time, end_time), dtype=int) - - event_trigger_data = [] - for trigger in trigger_idx: - target_idx = ind + trigger - - if np.min(target_idx) >= 0 and np.max(target_idx) < num_frame: - event_trigger_data.append(neural_data[target_idx]) - - # (num_event, cell_number, event_time_lambda) - event_trigger_data = np.array(event_trigger_data) - - return event_trigger_data - - -def ETA( - neural_data: FluoData, - behaviors_data: BehaviorData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(mean=TimeSeriesData): - function_id = output_dir.split("/")[-1] - print("start ETA:", function_id) - - neural_data = neural_data.data - behaviors_data = behaviors_data.data - - if params["transpose_x"]: - X = neural_data.transpose() - else: - X = neural_data - - if params["transpose_y"]: - Y = behaviors_data.transpose() - else: - Y = behaviors_data - - assert ( - X.shape[0] == Y.shape[0] - ), f""" - neural_data and behaviors_data is not same dimension, - neural.shape{X.shape}, behavior.shape{Y.shape}""" - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - cell_numbers = ind - X = X[:, ind] - - Y = Y[:, params["target_index"]] - - # calculate Triggers - trigger_idx = calc_trigger(Y, params["trigger_type"], params["trigger_threshold"]) - - # calculate Triggered average - event_trigger_data = calc_trigger_average( - X, trigger_idx, params["start_time"], params["end_time"] - ) - - # (cell_number, event_time_lambda) - if len(event_trigger_data) > 0: - mean = np.mean(event_trigger_data, axis=0) - sem = np.std(event_trigger_data, axis=0) / np.sqrt(len(event_trigger_data)) - mean = mean.transpose() - sem = sem.transpose() - else: - assert False, "Output data size is 0" - - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "mean": mean, - "sem": sem, - "num_sample": [len(mean)], +class ETAParams(BaseModel): + transpose_x: bool = Field(True) + transpose_y: bool = Field(False) + target_index: int = Field(1) + trigger_type: str = Field("up", description="available: up, down, cross") + trigger_threshold: int = Field(0) + start_time: int = Field(-10) + end_time: int = Field(10) + + +class ETA(AlgoTemplate): + def run( + self, + params: ETAParams, + neural_data: FluoData, + behaviors_data: BehaviorData, + iscell: IscellData = None, + ) -> dict(mean=TimeSeriesData): + print("start ETA:", self.function_id) + + neural_data = neural_data.data + behaviors_data = behaviors_data.data + + if params["transpose_x"]: + X = neural_data.transpose() + else: + X = neural_data + + if params["transpose_y"]: + Y = behaviors_data.transpose() + else: + Y = behaviors_data + + assert ( + X.shape[0] == Y.shape[0] + ), f""" + neural_data and behaviors_data is not same dimension, + neural.shape{X.shape}, behavior.shape{Y.shape}""" + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + cell_numbers = ind + X = X[:, ind] + + Y = Y[:, params["target_index"]] + + # calculate Triggers + trigger_idx = self.calc_trigger( + Y, params["trigger_type"], params["trigger_threshold"] + ) + + # calculate Triggered average + event_trigger_data = self.calc_trigger_average( + X, trigger_idx, params["start_time"], params["end_time"] + ) + + # (cell_number, event_time_lambda) + if len(event_trigger_data) > 0: + mean = np.mean(event_trigger_data, axis=0) + sem = np.std(event_trigger_data, axis=0) / np.sqrt(len(event_trigger_data)) + mean = mean.transpose() + sem = sem.transpose() + else: + assert False, "Output data size is 0" + + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "mean": mean, + "sem": sem, + "num_sample": [len(mean)], + } } - } - - min_value = np.min(mean, axis=1, keepdims=True) - max_value = np.max(mean, axis=1, keepdims=True) - norm_mean = (mean - min_value) / (max_value - min_value) - - info = {} - info["mean"] = TimeSeriesData( - mean, - std=sem, - index=list(np.arange(params["start_time"], params["end_time"])), - cell_numbers=cell_numbers if iscell is not None else None, - file_name="mean", - ) - info["mean_heatmap"] = HeatMapData( - norm_mean, - columns=list(np.arange(params["start_time"], params["end_time"])), - file_name="mean_heatmap", - ) - info["nwbfile"] = nwbfile - - return info + + min_value = np.min(mean, axis=1, keepdims=True) + max_value = np.max(mean, axis=1, keepdims=True) + norm_mean = (mean - min_value) / (max_value - min_value) + + info = {} + info["mean"] = TimeSeriesData( + mean, + std=sem, + index=list(np.arange(params["start_time"], params["end_time"])), + cell_numbers=cell_numbers if iscell is not None else None, + file_name="mean", + ) + info["mean_heatmap"] = HeatMapData( + norm_mean, + columns=list(np.arange(params["start_time"], params["end_time"])), + file_name="mean_heatmap", + ) + info["nwbfile"] = nwbfile + + return info + + @staticmethod + def calc_trigger(behavior_data, trigger_type, trigger_threshold): + flg = np.array(behavior_data > trigger_threshold, dtype=int) + if trigger_type == "up": + trigger_idx = np.where(np.ediff1d(flg) == 1)[0] + elif trigger_type == "down": + trigger_idx = np.where(np.ediff1d(flg) == -1)[0] + elif trigger_type == "cross": + trigger_idx = np.where(np.ediff1d(flg) != 0)[0] + else: + trigger_idx = np.where(np.ediff1d(flg) == 0)[0] + + return trigger_idx + + @staticmethod + def calc_trigger_average(neural_data, trigger_idx, start_time, end_time): + num_frame = neural_data.shape[0] + + ind = np.array(range(start_time, end_time), dtype=int) + + event_trigger_data = [] + for trigger in trigger_idx: + target_idx = ind + trigger + + if np.min(target_idx) >= 0 and np.max(target_idx) < num_frame: + event_trigger_data.append(neural_data[target_idx]) + + # (num_event, cell_number, event_time_lambda) + event_trigger_data = np.array(event_trigger_data) + + return event_trigger_data diff --git a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/params/cell_grouping.yaml b/studio/app/optinist/wrappers/optinist/basic_neural_analysis/params/cell_grouping.yaml deleted file mode 100644 index 295102c45..000000000 --- a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/params/cell_grouping.yaml +++ /dev/null @@ -1,4 +0,0 @@ -transpose: false -threshold: 1.0 -start_time: -10 -end_time: 0 diff --git a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/params/eta.yaml b/studio/app/optinist/wrappers/optinist/basic_neural_analysis/params/eta.yaml deleted file mode 100644 index 869a0f82e..000000000 --- a/studio/app/optinist/wrappers/optinist/basic_neural_analysis/params/eta.yaml +++ /dev/null @@ -1,9 +0,0 @@ -transpose_x: True -transpose_y: False - -target_index: 1 - -trigger_type: 'up' # ['up', 'down', 'cross'] -trigger_threshold: 0 -start_time: -10 -end_time: 10 diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/__init__.py b/studio/app/optinist/wrappers/optinist/dimension_reduction/__init__.py index 975f76778..17b182132 100644 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/__init__.py +++ b/studio/app/optinist/wrappers/optinist/dimension_reduction/__init__.py @@ -1,5 +1,5 @@ from studio.app.optinist.wrappers.optinist.dimension_reduction.cca import CCA -from studio.app.optinist.wrappers.optinist.dimension_reduction.dpca_fit import dpca_fit +from studio.app.optinist.wrappers.optinist.dimension_reduction.dpca_fit import DPCAFit from studio.app.optinist.wrappers.optinist.dimension_reduction.pca import PCA from studio.app.optinist.wrappers.optinist.dimension_reduction.tsne import TSNE @@ -9,7 +9,7 @@ "conda_name": "optinist", }, "dpca": { - "function": dpca_fit, + "function": DPCAFit, "conda_name": "optinist", }, "pca": { diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/cca.py b/studio/app/optinist/wrappers/optinist/dimension_reduction/cca.py index 8dab17355..54066fd3e 100644 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/cca.py +++ b/studio/app/optinist/wrappers/optinist/dimension_reduction/cca.py @@ -1,79 +1,101 @@ +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import BarData, ScatterData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import BehaviorData, FluoData, IscellData from studio.app.optinist.wrappers.optinist.utils import standard_norm -def CCA( - neural_data: FluoData, - behaviors_data: BehaviorData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - import numpy as np - from sklearn.cross_decomposition import CCA - - function_id = output_dir.split("/")[-1] - print("start cca:", function_id) - - neural_data = neural_data.data - behaviors_data = behaviors_data.data - - # data should be time x component matrix - if params["transpose_x"]: - X = neural_data.transpose() - else: - X = neural_data - - if params["transpose_y"]: - Y = behaviors_data.transpose() - else: - Y = behaviors_data - - assert ( - X.shape[0] == Y.shape[0] - ), f""" - neural_data and behaviors_data is not same dimension, - neural.shape{X.shape}, behavior.shape{Y.shape}""" - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[:, ind] - - Y = Y[:, params["target_index"]].reshape(-1, 1) - - # preprocessing - tX = standard_norm(X, params["standard_x_mean"], params["standard_x_std"]) - tY = standard_norm(Y, params["standard_y_mean"], params["standard_y_std"]) - - # calculate CCA - cca = CCA(**params["CCA"]) - projX, projY = cca.fit_transform(tX, tY) - - proj = np.concatenate([projX, projY], axis=1) - - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "projectedNd": proj, - "x_weights": cca.x_weights_, # singular vectors - "y_weights": cca.y_weights_, - "x_loadings_": cca.x_rotations_, - "y_loadings_": cca.x_rotations_, - "coef": cca.coef_, - "n_iter_": cca.n_iter_, - # 'n_features_in_': [cca.n_features_in_], +class CCAMainParams(BaseModel): + n_components: int = Field(2) + scale: bool = Field(True) + max_iter: int = Field(500) + tol: float = Field(0.000001) + copy_data: bool = Field(True, alias="copy") + + +class CCAParams(BaseModel): + standard_x_mean: bool = Field(True) + standard_x_std: bool = Field(True) + standard_y_mean: bool = Field(True) + standard_y_std: bool = Field(True) + transpose_x: bool = Field(True) + transpose_y: bool = Field(False) + target_index: int = Field(0) + CCA: CCAMainParams = Field(CCAMainParams()) + + +class CCA(AlgoTemplate): + def run( + self, + params: CCAParams, + neural_data: FluoData, + behaviors_data: BehaviorData, + iscell: IscellData = None, + ) -> dict(): + import numpy as np + from sklearn.cross_decomposition import CCA as SkCCA + + print("start cca:", self.function_id) + + neural_data = neural_data.data + behaviors_data = behaviors_data.data + + # data should be time x component matrix + if params["transpose_x"]: + X = neural_data.transpose() + else: + X = neural_data + + if params["transpose_y"]: + Y = behaviors_data.transpose() + else: + Y = behaviors_data + + assert ( + X.shape[0] == Y.shape[0] + ), f""" + neural_data and behaviors_data is not same dimension, + neural.shape{X.shape}, behavior.shape{Y.shape}""" + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[:, ind] + + Y = Y[:, params["target_index"]].reshape(-1, 1) + + # preprocessing + tX = standard_norm(X, params["standard_x_mean"], params["standard_x_std"]) + tY = standard_norm(Y, params["standard_y_mean"], params["standard_y_std"]) + + # calculate CCA + cca = SkCCA(**params["CCA"]) + projX, projY = cca.fit_transform(tX, tY) + + proj = np.concatenate([projX, projY], axis=1) + + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "projectedNd": proj, + "x_weights": cca.x_weights_, # singular vectors + "y_weights": cca.y_weights_, + "x_loadings_": cca.x_rotations_, + "y_loadings_": cca.x_rotations_, + "coef": cca.coef_, + "n_iter_": cca.n_iter_, + # 'n_features_in_': [cca.n_features_in_], + } } - } - info = { - "projectedNd": ScatterData(proj, file_name="projectedNd"), - "coef": BarData(cca.coef_.flatten(), file_name="coef"), - "nwbfile": nwbfile, - } + info = { + "projectedNd": ScatterData(proj, file_name="projectedNd"), + "coef": BarData(cca.coef_.flatten(), file_name="coef"), + "nwbfile": nwbfile, + } - return info + return info diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/dpca_fit.py b/studio/app/optinist/wrappers/optinist/dimension_reduction/dpca_fit.py index 1930e73a4..82cdd5396 100644 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/dpca_fit.py +++ b/studio/app/optinist/wrappers/optinist/dimension_reduction/dpca_fit.py @@ -1,242 +1,317 @@ +from typing import List, Optional, Union + import numpy as np import pandas as pd +from pydantic import BaseModel +from pydantic.dataclasses import Field +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import HeatMapData # , TimeSeriesData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import BehaviorData, FluoData, IscellData from studio.app.optinist.wrappers.optinist.utils import standard_norm -def calc_trigger(behavior_data, trigger_type, trigger_threshold): - # same function is also in the eta - flg = np.array(behavior_data > trigger_threshold, dtype=int) - if trigger_type == "up": - trigger_idx = np.where(np.ediff1d(flg) == 1)[0] - elif trigger_type == "down": - trigger_idx = np.where(np.ediff1d(flg) == -1)[0] - elif trigger_type == "cross": - trigger_idx = np.where(np.ediff1d(flg) != 0)[0] - else: - trigger_idx = np.where(np.ediff1d(flg) == 0)[0] - - return trigger_idx - - -def GetIndices(dims, outtype): - index = np.indices(dims) - ind = [] - for i in range(len(index)): - ind.append(index[i].flatten()) - ind = np.array(ind) - ind = ind.transpose() - - if outtype == "list": - out = [] - for i in range(ind.shape[1]): - out.append(ind[:, i]) - else: - out = ind - - return out - - -def createMatrix(D, triggers, stims, duration): - # D: num_timestamps x num_unit neural data - # triggers: index of trigger - # stims: list of stimulus property for each trigger - # durations: frames before and after trigger to use - - num_unit = D.shape[1] - num_triggers = len(triggers) - num_property = len(stims) - num_timepoints = duration[1] - duration[0] - - # X is the reshaped data for each trigger - X = np.zeros([num_triggers, num_unit, num_timepoints]) - for i in range(num_triggers): - X[i, :, :] = D[ - triggers[i] + duration[0] : triggers[i] + duration[1], : - ].transpose() - - # df is the table of stimulus conditions - # df: number of trigger x ( number of property +1 ) - # uq_stims = list of unique stimulus for each property - # num_uq_stims = number of unique_stims for each property - - columns = columns = list(map(str, range(num_property))) - columns.append("all") - df = pd.DataFrame(data=None, index=list(range(num_triggers)), columns=columns) - for i in range(num_property): - df[columns[i]] = stims[i] - - for i in range(num_triggers): - df.iloc[i, num_property] = "_".join(map(str, df.iloc[i, 0:2].values.tolist())) - - uq_list = df["all"].unique() - uq_stims = [] - num_uq_stims = [] - for i in range(num_property): - uq_stims.append(list(df.iloc[:, i].unique())) - num_uq_stims.append(len(uq_stims[i])) - - # check number of samples - # n: number of samples for each condition - # index: index of the trigger for each condition - # min_sample: minimum number of samples for a condition - n = np.zeros([len(uq_list)], dtype=int) - index = [] - for i in range(len(uq_list)): - index.append(list(df[df["all"] == uq_list[i]].index)) - n[i] = len(index[i]) - - min_sample = int(np.min(n)) - - # re-format the data (number of samples is set to the nun_sample) - X2 = np.zeros([min_sample, num_unit, num_timepoints] + num_uq_stims) - - for i in range(len(index)): - stims = list(df.iloc[index[i][0], 0 : df.shape[1] - 1]) - tgtind = [] - for j in range(len(stims)): - tgtind.append(uq_stims[j].index(stims[j])) - - if num_property == 1: - X2[0:min_sample, :, :, tgtind[0]] = X[index[i][0:min_sample], :, :] - elif num_property == 2: - X2[0:min_sample, :, :, tgtind[0], tgtind[1]] = X[ - index[i][0:min_sample], :, : - ] - elif num_property == 3: - X2[0:min_sample, :, :, tgtind[0], tgtind[1], tgtind[2]] = X[ - index[i][0:min_sample], :, : - ] - - else: - print("currently the number of condition category has to be less than 4") - return - - return X2 - - -def reshapeBehavior( - B, trigger_column, trigger_type, trigger_threshold, feature_columns -): - Trig = calc_trigger(B[:, trigger_column], trigger_type, trigger_threshold) - - features = [] - for i in range(len(feature_columns)): - features.append(B[Trig, feature_columns[i]]) - - return [Trig, features] - - -def dpca_fit( - neural_data: FluoData, - behaviors_data: BehaviorData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - # modules specific to function - from dPCA import dPCA - - function_id = output_dir.split("/")[-1] - print("start dpca:", function_id) - - neural_data = neural_data.data - behaviors_data = behaviors_data.data - - # neural data should be time x cells - if params["transpose"]: - X = neural_data.transpose() - else: - X = neural_data - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[:, ind] - - # preprocessing - X = standard_norm(X, params["standard_mean"], params["standard_std"]) - - # create trigger and features - [Trig, features] = reshapeBehavior( - behaviors_data, - params["trigger_column"], - params["trigger_type"], - params["trigger_threshold"], - params["feature_colums"], +class DPCAParams(BaseModel): + transpose: bool = Field(False) + standard_mean: bool = Field(True) + standard_std: bool = Field(True) + labels: Union[int, str] = Field( + "tbc", + description=( + "Labels of feature axis. " + "If int the corresponding number of labels are selected " + "from the alphabet 'abcde...'" + "--- seems int is not accepted. " + "numbers of characters should match the dimension of the data" + ), ) - X = createMatrix(X, Trig, features, params["trigger_duration"]) - - # calculate dPCA # - - # X: array - like, shape(n_samples, n_features_1, n_features_2, ...) - # Training data, where n_samples in the number of samples - # and n_features_j is the number - # of the j - features(where the axis correspond to different parameters). - - dpca = dPCA.dPCA( - labels=params["labels"], - join=params["join"], - regularizer=params["regularizer"], - n_components=params["n_components"], - copy=params["copy"], - n_iter=params["n_iter"], + join: Optional[dict] = Field( + None, + description=( + "Parameter combinations to join. " + "If a data set has parametrized by time t and stimulus s, " + "then dPCA will split the data into marginalizations " + "corresponding to 't', 's' and 'ts'. " + "At times, we want to join different marginalizations (like 's' and 'ts') " + "e.g. if we are only interested in the time-modulated stimulus components. " + "In this case, we would pass {'ts' : ['s','ts']}." + ), ) - - result = dpca.fit_transform(np.mean(X, axis=0), X) - keys = list(result.keys()) - - Out_forfigure = {} - # figure shows only assigned components and properties - for i in range(len(params["figure_features"])): - for j in range(len(params["figure_components"])): - tp = result[params["figure_features"][i]][ - params["figure_components"][j], - ] # 1st component - inds = GetIndices(tp.shape[1:], "matrix") - arr = np.zeros([tp.shape[0], inds.shape[0]]) - for m in range(inds.shape[0]): - for k in range(tp.shape[0]): - arr[k, m] = tp[tuple([k] + list(inds[m, :]))] - - Out_forfigure[ - params["figure_features"][i] - + "-component" - + str(params["figure_components"][j]) - ] = arr - Out_forfigure["features"] = list(Out_forfigure.keys()) - - names = [] - for i in range(inds.shape[0]): - names.append("feature" + "_".join(map(str, inds[i, :]))) - Out_forfigure["trace_names"] = names - - # NWB - tpdic = {} - for i in range(len(keys)): - tpdic[keys[i]] = result[keys[i]] - - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = {function_id: {**tpdic}} - - info = {} - for i in range(len(Out_forfigure["features"])): - # info[Out_forfigure["features"][i]] = TimeSeriesData( - # Out_forfigure[Out_forfigure["features"][i]].transpose(), - # std=None, - # index=None, - # file_name=Out_forfigure["features"][i], - # ) - - info[Out_forfigure["features"][i]] = HeatMapData( - Out_forfigure[Out_forfigure["features"][i]].transpose(), - columns=None, - file_name=Out_forfigure["features"][i], + regularizer: Union[float, str, None] = Field( + 0, + description=( + "Regularization parameter. This can take None, float, str 'auto'" + "If None or 0, then no regularization is applied. " + "For float, the regularization weight is regularizer*var(data). " + "If 'auto', the optimal regularization parameter is found " + "during fitting (might take some time)." + ), + ) + n_components: Union[int, dict, None] = Field( + 8, + description=( + "Number of components to keep." + "If n_components is int, then the same number of components are kept " + "in every marginalization. " + "Otherwise, the dict allows to set the number of components " + "in each marginalization (e.g. {'t' : 10, 'ts' : 5})." + ), + ) + copy_data: bool = Field( + True, + description=( + "If False, data passed to fit are overwritten and running fit(X). " + "transform(X) will not yield the expected results, " + "use fit_transform(X) instead." + ), + alias="copy", + ) + n_iter: int = Field( + 0, description="Number of iterations for randomized SVD solver (sklearn)." + ) + trigger_column: int = Field(1) + trigger_type: str = Field("down") + trigger_threshold: int = Field(0) + trigger_duration: List[int] = Field([-10, 10]) + feature_colums: List[int] = Field([3, 4]) + figure_components: List[int] = Field([0, 1]) + figure_features: List[str] = Field(["t", "b", "c", "tbc"]) + + +class DPCAFit(AlgoTemplate): + def run( + self, + params: DPCAParams, + neural_data: FluoData, + behaviors_data: BehaviorData, + iscell: IscellData = None, + ) -> dict(): + # modules specific to function + from dPCA import dPCA + + print("start dpca:", self.function_id) + + neural_data = neural_data.data + behaviors_data = behaviors_data.data + + # neural data should be time x cells + if params["transpose"]: + X = neural_data.transpose() + else: + X = neural_data + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[:, ind] + + # preprocessing + X = standard_norm(X, params["standard_mean"], params["standard_std"]) + + # create trigger and features + [Trig, features] = self.reshapeBehavior( + behaviors_data, + params["trigger_column"], + params["trigger_type"], + params["trigger_threshold"], + params["feature_colums"], + ) + X = self.createMatrix(X, Trig, features, params["trigger_duration"]) + + # calculate dPCA # + + # X: array - like, shape(n_samples, n_features_1, n_features_2, ...) + # Training data, where n_samples in the number of samples + # and n_features_j is the number + # of the j - features(where the axis correspond to different parameters). + + dpca = dPCA.dPCA( + labels=params["labels"], + join=params["join"], + regularizer=params["regularizer"], + n_components=params["n_components"], + copy=params["copy"], + n_iter=params["n_iter"], ) - info["nwbfile"] = nwbfile - return info + result = dpca.fit_transform(np.mean(X, axis=0), X) + keys = list(result.keys()) + + Out_forfigure = {} + # figure shows only assigned components and properties + for i in range(len(params["figure_features"])): + for j in range(len(params["figure_components"])): + tp = result[params["figure_features"][i]][ + params["figure_components"][j], + ] # 1st component + inds = self.GetIndices(tp.shape[1:], "matrix") + arr = np.zeros([tp.shape[0], inds.shape[0]]) + for m in range(inds.shape[0]): + for k in range(tp.shape[0]): + arr[k, m] = tp[tuple([k] + list(inds[m, :]))] + + Out_forfigure[ + params["figure_features"][i] + + "-component" + + str(params["figure_components"][j]) + ] = arr + Out_forfigure["features"] = list(Out_forfigure.keys()) + + names = [] + for i in range(inds.shape[0]): + names.append("feature" + "_".join(map(str, inds[i, :]))) + Out_forfigure["trace_names"] = names + + # NWB + tpdic = {} + for i in range(len(keys)): + tpdic[keys[i]] = result[keys[i]] + + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = {self.function_id: {**tpdic}} + + info = {} + for i in range(len(Out_forfigure["features"])): + # info[Out_forfigure["features"][i]] = TimeSeriesData( + # Out_forfigure[Out_forfigure["features"][i]].transpose(), + # std=None, + # index=None, + # file_name=Out_forfigure["features"][i], + # ) + + info[Out_forfigure["features"][i]] = HeatMapData( + Out_forfigure[Out_forfigure["features"][i]].transpose(), + columns=None, + file_name=Out_forfigure["features"][i], + ) + info["nwbfile"] = nwbfile + + return info + + @classmethod + def GetIndices(cls, dims, outtype): + index = np.indices(dims) + ind = [] + for i in range(len(index)): + ind.append(index[i].flatten()) + ind = np.array(ind) + ind = ind.transpose() + + if outtype == "list": + out = [] + for i in range(ind.shape[1]): + out.append(ind[:, i]) + else: + out = ind + + return out + + @classmethod + def reshapeBehavior( + cls, B, trigger_column, trigger_type, trigger_threshold, feature_columns + ): + Trig = cls.calc_trigger(B[:, trigger_column], trigger_type, trigger_threshold) + + features = [] + for i in range(len(feature_columns)): + features.append(B[Trig, feature_columns[i]]) + + return [Trig, features] + + @classmethod + def createMatrix(cls, D, triggers, stims, duration): + # D: num_timestamps x num_unit neural data + # triggers: index of trigger + # stims: list of stimulus property for each trigger + # durations: frames before and after trigger to use + + num_unit = D.shape[1] + num_triggers = len(triggers) + num_property = len(stims) + num_timepoints = duration[1] - duration[0] + + # X is the reshaped data for each trigger + X = np.zeros([num_triggers, num_unit, num_timepoints]) + for i in range(num_triggers): + X[i, :, :] = D[ + triggers[i] + duration[0] : triggers[i] + duration[1], : + ].transpose() + + # df is the table of stimulus conditions + # df: number of trigger x ( number of property +1 ) + # uq_stims = list of unique stimulus for each property + # num_uq_stims = number of unique_stims for each property + + columns = columns = list(map(str, range(num_property))) + columns.append("all") + df = pd.DataFrame(data=None, index=list(range(num_triggers)), columns=columns) + for i in range(num_property): + df[columns[i]] = stims[i] + + for i in range(num_triggers): + df.iloc[i, num_property] = "_".join( + map(str, df.iloc[i, 0:2].values.tolist()) + ) + + uq_list = df["all"].unique() + uq_stims = [] + num_uq_stims = [] + for i in range(num_property): + uq_stims.append(list(df.iloc[:, i].unique())) + num_uq_stims.append(len(uq_stims[i])) + + # check number of samples + # n: number of samples for each condition + # index: index of the trigger for each condition + # min_sample: minimum number of samples for a condition + n = np.zeros([len(uq_list)], dtype=int) + index = [] + for i in range(len(uq_list)): + index.append(list(df[df["all"] == uq_list[i]].index)) + n[i] = len(index[i]) + + min_sample = int(np.min(n)) + + # re-format the data (number of samples is set to the nun_sample) + X2 = np.zeros([min_sample, num_unit, num_timepoints] + num_uq_stims) + + for i in range(len(index)): + stims = list(df.iloc[index[i][0], 0 : df.shape[1] - 1]) + tgtind = [] + for j in range(len(stims)): + tgtind.append(uq_stims[j].index(stims[j])) + + if num_property == 1: + X2[0:min_sample, :, :, tgtind[0]] = X[index[i][0:min_sample], :, :] + elif num_property == 2: + X2[0:min_sample, :, :, tgtind[0], tgtind[1]] = X[ + index[i][0:min_sample], :, : + ] + elif num_property == 3: + X2[0:min_sample, :, :, tgtind[0], tgtind[1], tgtind[2]] = X[ + index[i][0:min_sample], :, : + ] + + else: + print( + "currently the number of condition category has to be less than 4" + ) + return + + return X2 + + @classmethod + def calc_trigger(cls, behavior_data, trigger_type, trigger_threshold): + # same function is also in the eta + flg = np.array(behavior_data > trigger_threshold, dtype=int) + if trigger_type == "up": + trigger_idx = np.where(np.ediff1d(flg) == 1)[0] + elif trigger_type == "down": + trigger_idx = np.where(np.ediff1d(flg) == -1)[0] + elif trigger_type == "cross": + trigger_idx = np.where(np.ediff1d(flg) != 0)[0] + else: + trigger_idx = np.where(np.ediff1d(flg) == 0)[0] + + return trigger_idx diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/params/cca.yaml b/studio/app/optinist/wrappers/optinist/dimension_reduction/params/cca.yaml deleted file mode 100644 index 022b84d24..000000000 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/params/cca.yaml +++ /dev/null @@ -1,17 +0,0 @@ -standard_x_mean: True -standard_x_std: True - -standard_y_mean: True -standard_y_std: True - -transpose_x: True -transpose_y: False - -target_index: 0 - -CCA: - n_components: 2 - scale: True - max_iter: 500 - tol: 0.000001 - copy: True diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/params/dpca.yaml b/studio/app/optinist/wrappers/optinist/dimension_reduction/params/dpca.yaml deleted file mode 100644 index 02f3dcd0b..000000000 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/params/dpca.yaml +++ /dev/null @@ -1,58 +0,0 @@ -transpose: False - -standard_mean: True -standard_std: True - - -labels: 'tbc' -# int or string -# Labels of feature axis. -# If int the corresponding number of labels are selected from the alphabet 'abcde...' -# --- seems int is not accepted. numbers of characters should match the dimension of the data -join: -# None or dict -# Parameter combinations to join -# If a data set has parametrized by time t and stimulus s, then dPCA will split -# the data into marginalizations corresponding to 't', 's' and 'ts'. At times, -# we want to join different marginalizations (like 's' and 'ts'), e.g. if -# we are only interested in the time-modulated stimulus components. In this case, -# we would pass {'ts' : ['s','ts']}. - -regularizer: 0 -# None, float, 'auto' -# Regularization parameter. If None or 0, then no regularization is applied. -# For float, the regularization weight is regularizer*var(data). If 'auto', the -# optimal regularization parameter is found during fitting (might take some time). - - -n_components: 8 -# None, int or dict -# Number of components to keep. -# If n_components is int, then the same number of components are kept in every -# marginalization. Otherwise, the dict allows to set the number of components -# in each marginalization (e.g. {'t' : 10, 'ts' : 5}). Defaults to 10. - -copy: True -# bool -# If False, data passed to fit are overwritten and running -# fit(X).transform(X) will not yield the expected results, -# use fit_transform(X) instead. - -n_iter: 0 -# int (default: 0) -# Number of iterations for randomized SVD solver (sklearn). - -trigger_column: 1 - -trigger_type: 'down' - -trigger_threshold: 0 - -trigger_duration: [-10, 10] - -feature_colums: [3,4] - - -figure_components: [0,1] - -figure_features: ['t', 'b', 'c', 'tbc'] diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/params/pca.yaml b/studio/app/optinist/wrappers/optinist/dimension_reduction/params/pca.yaml deleted file mode 100644 index 350a88b37..000000000 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/params/pca.yaml +++ /dev/null @@ -1,14 +0,0 @@ -# whether standardize the data or not -standard_mean: True -standard_std: True - -transpose: True - -PCA: - n_components: 2 - copy: True - whiten: False - svd_solver: 'auto' - tol: 0.0 - iterated_power: 'auto' - # random_state: 0 diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/params/tsne.yaml b/studio/app/optinist/wrappers/optinist/dimension_reduction/params/tsne.yaml deleted file mode 100644 index aa6eacc5a..000000000 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/params/tsne.yaml +++ /dev/null @@ -1,22 +0,0 @@ - -# whether standardize the data or not -standard_mean: True -standard_std: True - -transpose: True - -TSNE: - n_components: 2 - perplexity: 30.0 - early_exaggeration: 12.0 - learning_rate: 'warn' - n_iter: 1000 - n_iter_without_progress: 300 - min_grad_norm: 0.0000001 - metric: 'euclidean' - init: 'warn' - random_state: 0 - method: 'barnes_hut' - angle: 0.5 - n_jobs: 1 - square_distances: 'legacy' diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/pca.py b/studio/app/optinist/wrappers/optinist/dimension_reduction/pca.py index 2064a1ce1..bb7c4ac1d 100644 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/pca.py +++ b/studio/app/optinist/wrappers/optinist/dimension_reduction/pca.py @@ -1,75 +1,97 @@ +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import BarData, ScatterData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import FluoData, IscellData from studio.app.optinist.wrappers.optinist.utils import standard_norm -def PCA( - neural_data: FluoData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - # modules specific to function - import numpy as np - from sklearn.decomposition import PCA - - function_id = output_dir.split("/")[-1] - print("start PCA:", function_id) - - neural_data = neural_data.data - - # data should be time x component matrix - if params["transpose"]: - X = neural_data.transpose() - else: - X = neural_data - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[:, ind] - - # # preprocessing - tX = standard_norm(X, params["standard_mean"], params["standard_std"]) - - # calculate PCA - pca = PCA(**params["PCA"]) - proj_X = pca.fit_transform(tX) - - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "pca_projectedNd": proj_X, - "components": pca.components_, - "explained_variance": pca.explained_variance_, - "explained_variance_ratio": pca.explained_variance_ratio_, - "singular_values": pca.singular_values_, - "mean": pca.mean_, - "n_components": [pca.n_components_], - "n_samples": [pca.n_samples_], - "noise_variance": [pca.noise_variance_], - "n_features_in": [pca.n_features_in_], +class PCAMainParams(BaseModel): + n_components: int = Field(2) + copy_data: bool = Field(True, alias="copy") + whiten: bool = Field(False) + svd_solver: str = Field("auto") + tol: float = Field(0.0) + iterated_power: str = Field("auto") + # random_state: int = Field(0) + + +class PCAParams(BaseModel): + standard_mean: bool = Field(True) + standard_std: bool = Field(True) + transpose: bool = Field(True) + PCA: PCAMainParams = Field(PCAMainParams()) + + +class PCA(AlgoTemplate): + def run( + self, + params: PCAParams, + neural_data: FluoData, + iscell: IscellData = None, + ) -> dict(): + # modules specific to function + import numpy as np + from sklearn.decomposition import PCA + + print("start PCA:", self.function_id) + + neural_data = neural_data.data + + # data should be time x component matrix + if params["transpose"]: + X = neural_data.transpose() + else: + X = neural_data + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[:, ind] + + # # preprocessing + tX = standard_norm(X, params["standard_mean"], params["standard_std"]) + + # calculate PCA + pca = PCA(**params["PCA"]) + proj_X = pca.fit_transform(tX) + + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "pca_projectedNd": proj_X, + "components": pca.components_, + "explained_variance": pca.explained_variance_, + "explained_variance_ratio": pca.explained_variance_ratio_, + "singular_values": pca.singular_values_, + "mean": pca.mean_, + "n_components": [pca.n_components_], + "n_samples": [pca.n_samples_], + "noise_variance": [pca.noise_variance_], + "n_features_in": [pca.n_features_in_], + } } - } - - # import pdb; pdb.set_trace() - info = { - "explained_variance": BarData(pca.explained_variance_ratio_, file_name="evr"), - "projectedNd": ScatterData(proj_X, file_name="projectedNd"), - "contribution": BarData( - pca.components_, - index=[f"pca: {i}" for i in range(len(pca.components_))], - file_name="contribution", - ), - "cumsum_contribution": BarData( - np.cumsum(pca.components_, axis=0), - index=[f"pca: {i}" for i in range(len(pca.components_))], - file_name="cumsum_contribution", - ), - "nwbfile": nwbfile, - } - - return info + + # import pdb; pdb.set_trace() + info = { + "explained_variance": BarData( + pca.explained_variance_ratio_, file_name="evr" + ), + "projectedNd": ScatterData(proj_X, file_name="projectedNd"), + "contribution": BarData( + pca.components_, + index=[f"pca: {i}" for i in range(len(pca.components_))], + file_name="contribution", + ), + "cumsum_contribution": BarData( + np.cumsum(pca.components_, axis=0), + index=[f"pca: {i}" for i in range(len(pca.components_))], + file_name="cumsum_contribution", + ), + "nwbfile": nwbfile, + } + + return info diff --git a/studio/app/optinist/wrappers/optinist/dimension_reduction/tsne.py b/studio/app/optinist/wrappers/optinist/dimension_reduction/tsne.py index 402410103..81a5db7e8 100644 --- a/studio/app/optinist/wrappers/optinist/dimension_reduction/tsne.py +++ b/studio/app/optinist/wrappers/optinist/dimension_reduction/tsne.py @@ -1,50 +1,77 @@ +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import ScatterData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import FluoData, IscellData from studio.app.optinist.wrappers.optinist.utils import standard_norm -def TSNE( - neural_data: FluoData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - import numpy as np - from sklearn.manifold import TSNE +class TSNEMainParams(BaseModel): + n_components: int = Field(2) + perplexity: float = Field(30.0) + early_exaggeration: float = Field(12.0) + learning_rate: str = Field("warn") + n_iter: int = Field(1000) + n_iter_without_progress: int = Field(300) + min_grad_norm: float = Field(0.0000001) + metric: str = Field("euclidean") + init: str = Field("warn") + random_state: int = Field(0) + method: str = Field("barnes_hut") + angle: float = Field(0.5) + n_jobs: int = Field(1) + square_distances: str = Field("legacy") + + +class TSNEParams(BaseModel): + standard_mean: bool = Field(True) + standard_std: bool = Field(True) + transpose: bool = Field(True) + TSNE: TSNEMainParams = Field(TSNEMainParams()) + + +class TSNE(AlgoTemplate): + def run( + self, + params: TSNEParams, + neural_data: FluoData, + iscell: IscellData = None, + ) -> dict(): + import numpy as np + from sklearn.manifold import TSNE - function_id = output_dir.split("/")[-1] - print("start TSNE:", function_id) + print("start TSNE:", self.function_id) - neural_data = neural_data.data + neural_data = neural_data.data - # data should be time x component matrix - if params["transpose"]: - X = neural_data.transpose() - else: - X = neural_data + # data should be time x component matrix + if params["transpose"]: + X = neural_data.transpose() + else: + X = neural_data - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[:, ind] + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[:, ind] - # preprocessing - tX = standard_norm(X, params["standard_mean"], params["standard_std"]) + # preprocessing + tX = standard_norm(X, params["standard_mean"], params["standard_std"]) - # calculate TSNE - tsne = TSNE(**params["TSNE"]) + # calculate TSNE + tsne = TSNE(**params["TSNE"]) - proj_X = tsne.fit_transform(tX) + proj_X = tsne.fit_transform(tX) - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = {function_id: {"projectedNd": proj_X}} + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = {self.function_id: {"projectedNd": proj_X}} - info = { - "projectedNd": ScatterData(proj_X, file_name="projectedNd"), - "nwbfile": nwbfile, - } + info = { + "projectedNd": ScatterData(proj_X, file_name="projectedNd"), + "nwbfile": nwbfile, + } - return info + return info diff --git a/studio/app/optinist/wrappers/optinist/neural_decoding/glm.py b/studio/app/optinist/wrappers/optinist/neural_decoding/glm.py index ca6f87d7c..48c1f84ab 100644 --- a/studio/app/optinist/wrappers/optinist/neural_decoding/glm.py +++ b/studio/app/optinist/wrappers/optinist/neural_decoding/glm.py @@ -4,104 +4,129 @@ # # https://www.statsmodels.org/stable/glm.html +from typing import List, Optional + +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import BarData, HTMLData, ScatterData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import BehaviorData, FluoData, IscellData from studio.app.optinist.wrappers.optinist.utils import standard_norm -def GLM( - neural_data: FluoData, - behaviors_data: BehaviorData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - # modules specific to function - import numpy as np - import pandas as pd - import statsmodels.api as sm - - function_id = output_dir.split("/")[-1] - print("start glm:", function_id) - - neural_data = neural_data.data - behaviors_data = behaviors_data.data - - # data should be time x component matrix - if params["transpose_x"]: - X = neural_data.transpose() - else: - X = neural_data - - if params["transpose_y"]: - Y = behaviors_data.transpose() - else: - Y = behaviors_data - - assert ( - X.shape[0] == Y.shape[0] - ), f""" - neural_data and behaviors_data is not same dimension, - neural.shape{X.shape}, behavior.shape{Y.shape}""" - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[:, ind] - Y = Y[:, ind] - - Y = Y[:, params["target_index"]].reshape(-1, 1) - - # preprocessing - tX = standard_norm(X, params["standard_x_mean"], params["standard_x_std"]) - tY = standard_norm(Y, params["standard_y_mean"], params["standard_y_std"]) - - # calculate GLM - tX = pd.DataFrame(tX) - tY = pd.DataFrame(tY) - - if params["add_constant"]: - tX = sm.add_constant(tX, prepend=False) - - # set family - link = getattr(sm.genmod.families.links, params["link"])() # noqa: F841 - family = eval(f"sm.families.{params['family']}(link=link)") - - # model fit - model = sm.GLM(tY, tX, family=family, **params["GLM"]) - Res = model.fit() - - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "actual_predicted": np.array([Res._endog, Res.mu]).transpose(), - "params": Res.params.values, - "pvalues": Res.pvalues.values, - "tvalues": Res.tvalues.values, # z - "aic": [Res.aic], - "bic_llf": [Res.bic_llf], - "llf": [Res.llf], # log-Likelihood - "pearson_chi2": [Res.pearson_chi2], - "df_model": [Res.df_model], - "df_resid": [Res.df_resid], - "scale": [Res.scale], - "mu": Res.mu, - "endog": Res._endog, +class GLMMainParams(BaseModel): + offset: Optional[List] = Field(None) + exposure: Optional[List] = Field(None) + missing: Optional[str] = Field(None) + + +class GLMParams(BaseModel): + standard_x_mean: bool = Field(True) + standard_x_std: bool = Field(True) + standard_y_mean: bool = Field(True) + standard_y_std: bool = Field(True) + transpose_x: bool = Field(True) + transpose_y: bool = Field(False) + target_index: int = Field(0) + add_constant: bool = Field(False) + link: str = Field("log") + family: str = Field("Gaussian") + GLM: GLMMainParams = Field(GLMMainParams()) + + +class GLM(AlgoTemplate): + def run( + self, + params: GLMParams, + neural_data: FluoData, + behaviors_data: BehaviorData, + iscell: IscellData = None, + ) -> dict(): + # modules specific to function + import numpy as np + import pandas as pd + import statsmodels.api as sm + + print("start glm:", self.function_id) + + neural_data = neural_data.data + behaviors_data = behaviors_data.data + + # data should be time x component matrix + if params["transpose_x"]: + X = neural_data.transpose() + else: + X = neural_data + + if params["transpose_y"]: + Y = behaviors_data.transpose() + else: + Y = behaviors_data + + assert ( + X.shape[0] == Y.shape[0] + ), f""" + neural_data and behaviors_data is not same dimension, + neural.shape{X.shape}, behavior.shape{Y.shape}""" + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[:, ind] + Y = Y[:, ind] + + Y = Y[:, params["target_index"]].reshape(-1, 1) + + # preprocessing + tX = standard_norm(X, params["standard_x_mean"], params["standard_x_std"]) + tY = standard_norm(Y, params["standard_y_mean"], params["standard_y_std"]) + + # calculate GLM + tX = pd.DataFrame(tX) + tY = pd.DataFrame(tY) + + if params["add_constant"]: + tX = sm.add_constant(tX, prepend=False) + + # set family + link = getattr(sm.genmod.families.links, params["link"])() # noqa: F841 + family = eval(f"sm.families.{params['family']}(link=link)") + + # model fit + model = sm.GLM(tY, tX, family=family, **params["GLM"]) + Res = model.fit() + + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "actual_predicted": np.array([Res._endog, Res.mu]).transpose(), + "params": Res.params.values, + "pvalues": Res.pvalues.values, + "tvalues": Res.tvalues.values, # z + "aic": [Res.aic], + "bic_llf": [Res.bic_llf], + "llf": [Res.llf], # log-Likelihood + "pearson_chi2": [Res.pearson_chi2], + "df_model": [Res.df_model], + "df_resid": [Res.df_resid], + "scale": [Res.scale], + "mu": Res.mu, + "endog": Res._endog, + } } - } - - # main results for plot - # plot should be reconsidered --- what they should be! - info = { - "actual_predicted": ScatterData( - np.array([Res._endog, Res.mu]).transpose(), file_name="actual_predicted" - ), - "params": BarData(Res.params.values, file_name="params"), - "textout": HTMLData(Res.summary().as_html(), file_name="textout"), - "nwbfile": nwbfile, - } - - return info + + # main results for plot + # plot should be reconsidered --- what they should be! + info = { + "actual_predicted": ScatterData( + np.array([Res._endog, Res.mu]).transpose(), file_name="actual_predicted" + ), + "params": BarData(Res.params.values, file_name="params"), + "textout": HTMLData(Res.summary().as_html(), file_name="textout"), + "nwbfile": nwbfile, + } + + return info diff --git a/studio/app/optinist/wrappers/optinist/neural_decoding/lda.py b/studio/app/optinist/wrappers/optinist/neural_decoding/lda.py index 5a3a8bc96..cb06cd1f3 100644 --- a/studio/app/optinist/wrappers/optinist/neural_decoding/lda.py +++ b/studio/app/optinist/wrappers/optinist/neural_decoding/lda.py @@ -1,80 +1,111 @@ +from typing import List, Optional, Union + +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import BarData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import BehaviorData, FluoData, IscellData from studio.app.optinist.wrappers.optinist.utils import standard_norm -def LDA( - neural_data: FluoData, - behaviors_data: BehaviorData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - # modules specific to function - import numpy as np - from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA - from sklearn.model_selection import StratifiedKFold - - function_id = output_dir.split("/")[-1] - print("start LDA:", function_id) - - neural_data = neural_data.data - behaviors_data = behaviors_data.data - - # data should be time x component matrix - if params["transpose_x"]: - X = neural_data.transpose() - else: - X = neural_data - - if params["transpose_y"]: - Y = behaviors_data.transpose() - else: - Y = behaviors_data - - assert ( - X.shape[0] == Y.shape[0] - ), f""" - neural_data and behaviors_data is not same dimension, - neural.shape{neural_data.shape}, behavior.shape{behaviors_data.shape}""" - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[:, ind] - - Y = Y[:, params["target_index"]].reshape(-1, 1) - - # preprocessing - tX = standard_norm(X, params["standard_x_mean"], params["standard_x_std"]) - - # cross validation of LDA model - skf = StratifiedKFold(**params["CV"]) - - score = [] - classifier = [] - for train_index, test_index in skf.split(tX, Y): - clf = LDA(**params["LDA"]) - - if tX.shape[0] == 1: - clf.fit(tX[train_index].reshape(-1, 1), Y[train_index]) - score.append(clf.score(tX[test_index].reshape(-1, 1), Y[test_index])) - classifier.append(clf) +class LDACVParams(BaseModel): + n_splits: int = Field(5, description="Number of folds. Must be at least 2.") + shuffle: bool = Field(False) + # random_state: int = Field(0) + + +class LDAMainParams(BaseModel): + solver: str = Field("svd") + shrinkage: Union[str, float, None] = Field(None) + priors: Optional[List] = Field(None) + n_components: Optional[int] = Field(None) + store_covariance: bool = Field(False) + tol: float = Field(0.0001) + covariance_estimator: Optional[str] = Field(None) + + +class LDAParams(BaseModel): + standard_x_mean: bool = Field(True) + standard_x_std: bool = Field(True) + transpose_x: bool = Field(True) + transpose_y: bool = Field(False) + target_index: int = Field(1) + CV: LDACVParams = Field(LDACVParams()) + LDA: LDAMainParams = Field(LDAMainParams()) + + +class LDA(AlgoTemplate): + def run( + self, + params: LDAParams, + neural_data: FluoData, + behaviors_data: BehaviorData, + iscell: IscellData = None, + ) -> dict(): + # modules specific to function + import numpy as np + from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA + from sklearn.model_selection import StratifiedKFold + + print("start LDA:", self.function_id) + + neural_data = neural_data.data + behaviors_data = behaviors_data.data + + # data should be time x component matrix + if params["transpose_x"]: + X = neural_data.transpose() + else: + X = neural_data + + if params["transpose_y"]: + Y = behaviors_data.transpose() else: - clf.fit(tX[train_index, :], Y[train_index]) - score.append(clf.score(tX[test_index, :], Y[test_index])) - classifier.append(clf) - - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "score": score, + Y = behaviors_data + + assert ( + X.shape[0] == Y.shape[0] + ), f""" + neural_data and behaviors_data is not same dimension, + neural.shape{neural_data.shape}, behavior.shape{behaviors_data.shape}""" + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[:, ind] + + Y = Y[:, params["target_index"]].reshape(-1, 1) + + # preprocessing + tX = standard_norm(X, params["standard_x_mean"], params["standard_x_std"]) + + # cross validation of LDA model + skf = StratifiedKFold(**params["CV"]) + + score = [] + classifier = [] + for train_index, test_index in skf.split(tX, Y): + clf = LDA(**params["LDA"]) + + if tX.shape[0] == 1: + clf.fit(tX[train_index].reshape(-1, 1), Y[train_index]) + score.append(clf.score(tX[test_index].reshape(-1, 1), Y[test_index])) + classifier.append(clf) + else: + clf.fit(tX[train_index, :], Y[train_index]) + score.append(clf.score(tX[test_index, :], Y[test_index])) + classifier.append(clf) + + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "score": score, + } } - } - info = {"score": BarData(score, file_name="score"), "nwbfile": nwbfile} + info = {"score": BarData(score, file_name="score"), "nwbfile": nwbfile} - return info + return info diff --git a/studio/app/optinist/wrappers/optinist/neural_decoding/params/glm.yaml b/studio/app/optinist/wrappers/optinist/neural_decoding/params/glm.yaml deleted file mode 100644 index 515d7595c..000000000 --- a/studio/app/optinist/wrappers/optinist/neural_decoding/params/glm.yaml +++ /dev/null @@ -1,19 +0,0 @@ -standard_x_mean: True -standard_x_std: True - -standard_y_mean: True -standard_y_std: True - -transpose_x: True -transpose_y: False - -target_index: 0 - -add_constant: False -link: 'log' -family: 'Gaussian' - -GLM: - offset: - exposure: - missing: diff --git a/studio/app/optinist/wrappers/optinist/neural_decoding/params/lda.yaml b/studio/app/optinist/wrappers/optinist/neural_decoding/params/lda.yaml deleted file mode 100644 index fe4540a5a..000000000 --- a/studio/app/optinist/wrappers/optinist/neural_decoding/params/lda.yaml +++ /dev/null @@ -1,31 +0,0 @@ -# whether standardize the data or not -standard_x_mean: True -standard_x_std: True - -transpose_x: True -transpose_y: False - -target_index: 1 - -############################## -# StratifiedKFold parameters -############################## -# n_splits = int, default=5 -# Number of folds. Must be at least 2. -CV: - n_splits: 5 - shuffle: False - # random_state: 0 - - -############################## -# LDA parameters -############################## -LDA: - solver: 'svd' - shrinkage: - priors: - n_components: null - store_covariance: False - tol: 0.0001 - covariance_estimator: diff --git a/studio/app/optinist/wrappers/optinist/neural_decoding/params/svm.yaml b/studio/app/optinist/wrappers/optinist/neural_decoding/params/svm.yaml deleted file mode 100644 index 27668ff86..000000000 --- a/studio/app/optinist/wrappers/optinist/neural_decoding/params/svm.yaml +++ /dev/null @@ -1,70 +0,0 @@ - - -# whether standardize the data or not -standard_x_mean: True -standard_x_std: True - -transpose_x: True -transpose_y: False - -target_index: 1 - -############################## -# grid search parameters -############################## - -# whether grid search is performed or not -use_grid_search: True - -# grids for grid search -grid_search: - param_grid: - C: [0.001, 0.01, 0.1] - kernel: ['linear'] - degree: [3] - gamma: ['scale'] - coef0: [0.0] - shrinking: [True] - tol: [0.001] - decision_function_shape: ['ovr'] - - CV: - scoring: 'accuracy' - n_jobs: 1 - refit: True - cv: - verbose: 3 - pre_dispatch: '2*n_jobs' - error_score: -1000000 - return_train_score: False - - -############################## -# StratifiedKFold parameters -############################## -CV: - n_splits: 5 - shuffle: True - # random_state: 0 - -############################## -# SVC parameters -############################## - -# C= float, default=1.0 -# Regularization parameter. The strength of the regularization is inversely proportional to C. Must be strictly positive. The penalty is a squared l2 penalty. -SVC: - C: 1.0 - kernel: 'rbf' - degree: 3 - gamma: 'scale' - coef0: 0.0 - shrinking: True - probability: False - tol: 0.001 - cache_size: 200 - class_weight: - max_iter: -1 - decision_function_shape: 'ovr' - break_ties: False - random_state: 0 diff --git a/studio/app/optinist/wrappers/optinist/neural_decoding/svm.py b/studio/app/optinist/wrappers/optinist/neural_decoding/svm.py index 4fad22105..ddb86f009 100644 --- a/studio/app/optinist/wrappers/optinist/neural_decoding/svm.py +++ b/studio/app/optinist/wrappers/optinist/neural_decoding/svm.py @@ -1,97 +1,164 @@ +from typing import List, Optional, Union + +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import BarData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import BehaviorData, FluoData, IscellData from studio.app.optinist.wrappers.optinist.utils import standard_norm -def SVM( - neural_data: FluoData, - behaviors_data: BehaviorData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - # modules specific to function - import numpy as np - from sklearn import svm - from sklearn.model_selection import GridSearchCV, StratifiedKFold - - function_id = output_dir.split("/")[-1] - print("start SVM:", function_id) - - neural_data = neural_data.data - behaviors_data = behaviors_data.data - - # data should be time x component matrix - if params["transpose_x"]: - X = neural_data.transpose() - else: - X = neural_data - - if params["transpose_y"]: - Y = behaviors_data.transpose() - else: - Y = behaviors_data - - assert ( - X.shape[0] == Y.shape[0] - ), f""" - neural_data and behaviors_data is not same dimension, - neural.shape{X.shape}, behavior.shape{Y.shape}""" - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[:, ind] - - Y = Y[:, params["target_index"]].reshape(-1, 1) - - # preprocessing - tX = standard_norm(X, params["standard_x_mean"], params["standard_x_std"]) - - hp = params["SVC"] - - # SVM determination of hyper parameters if needed - gs_clf = [] - if params["use_grid_search"]: - param_grid = [params["grid_search"]["param_grid"]] - gs_clf = GridSearchCV( - svm.SVC(), param_grid=param_grid, **params["grid_search"]["CV"] - ) - - gs_clf.fit(tX, Y) - - # insert best parameter to hp dictionary - keys = list(gs_clf.best_params_.keys()) - for i in range(len(keys)): - hp[keys[i]] = gs_clf.best_params_[keys[i]] - - # cross validation of SVM using best grid search paraneters - skf = StratifiedKFold(**params["CV"]) - - score = [] - classifier = [] - for train_index, test_index in skf.split(tX, Y): - clf = svm.SVC(**hp) - - if tX.shape[0] == 1: - clf.fit(tX[train_index].reshape(-1, 1), Y[train_index]) - score.append(clf.score(tX[test_index].reshape(-1, 1), Y[test_index])) +class SVMGridSearchParamGridParams(BaseModel): + C: List[float] = Field([0.001, 0.01, 0.1]) + kernel: List[str] = Field(["linear"]) + degree: List[int] = Field([3]) + gamma: List[str] = Field(["scale"]) + coef0: List[float] = Field([0.0]) + shrinking: List[bool] = Field([True]) + tol: List[float] = Field([0.001]) + decision_function_shape: List[str] = Field(["ovr"]) + + +class SVMGridSearchCVParams(BaseModel): + scoring: str = Field("accuracy") + n_jobs: int = Field(1) + refit: bool = Field(True) + cv: Optional[int] = Field(None) + verbose: int = Field(3) + pre_dispatch: str = Field("2*n_jobs") + error_score: int = Field(-1000000) + return_train_score: bool = Field(False) + + +class SVMGridSearchParams(BaseModel): + param_grid: SVMGridSearchParamGridParams = Field(SVMGridSearchParamGridParams()) + CV: SVMGridSearchCVParams = Field(SVMGridSearchCVParams()) + + +class SVMCVParams(BaseModel): + n_splits: int = Field(5) + shuffle: bool = Field(True) + # random_state: int = Field(0) + + +class SVMSVCParams(BaseModel): + C: float = Field(1.0) + kernel: str = Field("rbf") + degree: int = Field(3) + gamma: str = Field("scale") + coef0: float = Field(0.0) + shrinking: bool = Field(True) + probability: bool = Field(False) + tol: float = Field(0.001) + cache_size: float = Field(200) + class_weight: Union[str, dict, None] = Field(None) + max_iter: int = Field(-1) + decision_function_shape: str = Field("ovr") + break_ties: bool = Field(False) + random_state: int = Field(0) + + +class SVMParams(BaseModel): + standard_x_mean: bool = Field(True) + standard_x_std: bool = Field(True) + transpose_x: bool = Field(True) + transpose_y: bool = Field(False) + target_index: int = Field(1) + use_grid_search: bool = Field(True) + grid_search: SVMGridSearchParams = Field(SVMGridSearchParams()) + CV: SVMCVParams = Field(SVMCVParams()) + SVC: SVMSVCParams = Field(SVMSVCParams()) + + +class SVM(AlgoTemplate): + def run( + self, + params: SVMParams, + neural_data: FluoData, + behaviors_data: BehaviorData, + iscell: IscellData = None, + ) -> dict(): + # modules specific to function + import numpy as np + from sklearn import svm + from sklearn.model_selection import GridSearchCV, StratifiedKFold + + print("start SVM:", self.function_id) + + neural_data = neural_data.data + behaviors_data = behaviors_data.data + + # data should be time x component matrix + if params["transpose_x"]: + X = neural_data.transpose() else: - clf.fit(tX[train_index, :], Y[train_index]) - score.append(clf.score(tX[test_index, :], Y[test_index])) + X = neural_data - classifier.append(clf) - - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "score": score, + if params["transpose_y"]: + Y = behaviors_data.transpose() + else: + Y = behaviors_data + + assert ( + X.shape[0] == Y.shape[0] + ), f""" + neural_data and behaviors_data is not same dimension, + neural.shape{X.shape}, behavior.shape{Y.shape}""" + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[:, ind] + + Y = Y[:, params["target_index"]].reshape(-1, 1) + + # preprocessing + tX = standard_norm(X, params["standard_x_mean"], params["standard_x_std"]) + + hp = params["SVC"] + + # SVM determination of hyper parameters if needed + gs_clf = [] + if params["use_grid_search"]: + param_grid = [params["grid_search"]["param_grid"]] + gs_clf = GridSearchCV( + svm.SVC(), param_grid=param_grid, **params["grid_search"]["CV"] + ) + + gs_clf.fit(tX, Y) + + # insert best parameter to hp dictionary + keys = list(gs_clf.best_params_.keys()) + for i in range(len(keys)): + hp[keys[i]] = gs_clf.best_params_[keys[i]] + + # cross validation of SVM using best grid search paraneters + skf = StratifiedKFold(**params["CV"]) + + score = [] + classifier = [] + for train_index, test_index in skf.split(tX, Y): + clf = svm.SVC(**hp) + + if tX.shape[0] == 1: + clf.fit(tX[train_index].reshape(-1, 1), Y[train_index]) + score.append(clf.score(tX[test_index].reshape(-1, 1), Y[test_index])) + else: + clf.fit(tX[train_index, :], Y[train_index]) + score.append(clf.score(tX[test_index, :], Y[test_index])) + + classifier.append(clf) + + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "score": score, + } } - } - info = {"score": BarData(score, file_name="score"), "nwbfile": nwbfile} + info = {"score": BarData(score, file_name="score"), "nwbfile": nwbfile} - return info + return info diff --git a/studio/app/optinist/wrappers/optinist/neural_population_analysis/__init__.py b/studio/app/optinist/wrappers/optinist/neural_population_analysis/__init__.py index f0772db35..6aec58bb2 100644 --- a/studio/app/optinist/wrappers/optinist/neural_population_analysis/__init__.py +++ b/studio/app/optinist/wrappers/optinist/neural_population_analysis/__init__.py @@ -1,8 +1,8 @@ from studio.app.optinist.wrappers.optinist.neural_population_analysis.correlation import ( # noqa: E501 - correlation, + Correlation, ) from studio.app.optinist.wrappers.optinist.neural_population_analysis.cross_correlation import ( # noqa: E501 - cross_correlation, + CrossCorrelation, ) from studio.app.optinist.wrappers.optinist.neural_population_analysis.granger import ( Granger, @@ -10,10 +10,10 @@ neural_population_analysis_wrapper_dict = { "correlation": { - "function": correlation, + "function": Correlation, }, "cross_correlation": { - "function": cross_correlation, + "function": CrossCorrelation, "conda_name": "optinist", }, "granger": { diff --git a/studio/app/optinist/wrappers/optinist/neural_population_analysis/correlation.py b/studio/app/optinist/wrappers/optinist/neural_population_analysis/correlation.py index e9e4c91ab..9b595fb4a 100644 --- a/studio/app/optinist/wrappers/optinist/neural_population_analysis/correlation.py +++ b/studio/app/optinist/wrappers/optinist/neural_population_analysis/correlation.py @@ -1,51 +1,58 @@ +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import HeatMapData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import FluoData, IscellData -def correlation( - neural_data: FluoData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - import numpy as np +class CorrelationParams(BaseModel): + transpose: bool = Field(True) + - function_id = output_dir.split("/")[-1] - print("start correlation:", function_id) +class Correlation(AlgoTemplate): + def run( + self, + params: CorrelationParams, + neural_data: FluoData, + iscell: IscellData = None, + ) -> dict(): + import numpy as np - neural_data = neural_data.data + print("start correlation:", self.function_id) - # data should be time x component matrix - if params["transpose"]: - X = neural_data.transpose() - else: - X = neural_data + neural_data = neural_data.data - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[ind, :] + # data should be time x component matrix + if params["transpose"]: + X = neural_data.transpose() + else: + X = neural_data - num_cell = X.shape[0] + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[ind, :] - # calculate correlation - corr = np.corrcoef(X) - for i in range(num_cell): - corr[i, i] = np.nan + num_cell = X.shape[0] - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "corr": corr, + # calculate correlation + corr = np.corrcoef(X) + for i in range(num_cell): + corr[i, i] = np.nan + + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "corr": corr, + } } - } - info = { - "corr": HeatMapData(corr, file_name="corr"), - "nwbfile": nwbfile, - } + info = { + "corr": HeatMapData(corr, file_name="corr"), + "nwbfile": nwbfile, + } - return info + return info diff --git a/studio/app/optinist/wrappers/optinist/neural_population_analysis/cross_correlation.py b/studio/app/optinist/wrappers/optinist/neural_population_analysis/cross_correlation.py index a188b6293..d73a046c2 100644 --- a/studio/app/optinist/wrappers/optinist/neural_population_analysis/cross_correlation.py +++ b/studio/app/optinist/wrappers/optinist/neural_population_analysis/cross_correlation.py @@ -1,116 +1,134 @@ +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import TimeSeriesData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import FluoData, IscellData -def cross_correlation( - neural_data: FluoData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - import itertools - - import numpy as np - import scipy.signal as ss - import scipy.stats as stats - from tqdm import tqdm - - function_id = output_dir.split("/")[-1] - print("start cross_correlation:", function_id) - - neural_data = neural_data.data - - # data should be time x component matrix - if params["transpose"]: - X = neural_data.transpose() - else: - X = neural_data - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[ind, :] - - # calculate cross correlation - num_cell = X.shape[0] - data_len = X.shape[1] - shuffle_num = params["shuffle_sample_number"] - - lags = ss.correlation_lags(data_len, data_len, mode="same") - ind = np.intersect1d( - np.where(-params["lags"] <= lags), np.where(lags <= params["lags"]) +class CrossCorrelationParams(BaseModel): + transpose: bool = Field(True) + method: str = Field("direct") + lags: int = Field(100, description="number of frames to show (+-frames)") + shuffle_sample_number: int = Field( + 100, description="number of cross correlation caliculation for randomization" ) - x = lags[ind] + shuffle_confidence_interval: float = Field(0.95) + + +class CrossCorrelation(AlgoTemplate): + def run( + self, + params: CrossCorrelationParams, + neural_data: FluoData, + iscell: IscellData = None, + ) -> dict(): + import itertools + + import numpy as np + import scipy.signal as ss + import scipy.stats as stats + from tqdm import tqdm + + print("start cross_correlation:", self.function_id) + + neural_data = neural_data.data + + # data should be time x component matrix + if params["transpose"]: + X = neural_data.transpose() + else: + X = neural_data + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[ind, :] + + # calculate cross correlation + num_cell = X.shape[0] + data_len = X.shape[1] + shuffle_num = params["shuffle_sample_number"] + + lags = ss.correlation_lags(data_len, data_len, mode="same") + ind = np.intersect1d( + np.where(-params["lags"] <= lags), np.where(lags <= params["lags"]) + ) + x = lags[ind] - mat = np.zeros([num_cell, num_cell, len(x)]) - s_confint = np.zeros([num_cell, num_cell, len(x), 2]) - s_mean = np.zeros([num_cell, num_cell, len(x)]) + mat = np.zeros([num_cell, num_cell, len(x)]) + s_confint = np.zeros([num_cell, num_cell, len(x), 2]) + s_mean = np.zeros([num_cell, num_cell, len(x)]) - for i in tqdm(range(num_cell)): - # for i in tqdm(range(3)): - for j in tqdm(range(num_cell)): - ccvals = ss.correlate( - X[i, :], X[j, :], method=params["method"], mode="same" - ) - mat[i, j, :] = ccvals[ind] - - # baseline - tp_s_mat = np.zeros([len(x), shuffle_num]) - for k in range(shuffle_num): - tp = X[j, :] - np.random.shuffle(tp) - ccvals = ss.correlate(X[i, :], tp, method=params["method"], mode="same") - tp_s_mat[:, k] = ccvals[ind] - - b_mean = np.mean(tp_s_mat, axis=1) - b_sem = stats.sem(tp_s_mat, axis=1) - b_df = shuffle_num - 1 - - interval = np.zeros([len(x), 2]) - for k in range(len(x)): - interval[k, :] = stats.t.interval( - params["shuffle_confidence_interval"], b_df, loc=0, scale=b_sem[k] + for i in tqdm(range(num_cell)): + # for i in tqdm(range(3)): + for j in tqdm(range(num_cell)): + ccvals = ss.correlate( + X[i, :], X[j, :], method=params["method"], mode="same" ) + mat[i, j, :] = ccvals[ind] + + # baseline + tp_s_mat = np.zeros([len(x), shuffle_num]) + for k in range(shuffle_num): + tp = X[j, :] + np.random.shuffle(tp) + ccvals = ss.correlate( + X[i, :], tp, method=params["method"], mode="same" + ) + tp_s_mat[:, k] = ccvals[ind] + + b_mean = np.mean(tp_s_mat, axis=1) + b_sem = stats.sem(tp_s_mat, axis=1) + b_df = shuffle_num - 1 + + interval = np.zeros([len(x), 2]) + for k in range(len(x)): + interval[k, :] = stats.t.interval( + params["shuffle_confidence_interval"], + b_df, + loc=0, + scale=b_sem[k], + ) + + s_confint[i, j, :, 0] = interval[:, 0] + s_confint[i, j, :, 1] = interval[:, 1] + s_mean[i, j, :] = b_mean + + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "mat": mat, + "baseline": s_mean, + "base_confint": s_confint, + } + } - s_confint[i, j, :, 0] = interval[:, 0] - s_confint[i, j, :, 1] = interval[:, 1] - s_mean[i, j, :] = b_mean - - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "mat": mat, - "baseline": s_mean, - "base_confint": s_confint, + info = { + # 'cross_correlation': ScatterData(mat), + "nwbfile": nwbfile } - } - - info = { - # 'cross_correlation': ScatterData(mat), - "nwbfile": nwbfile - } - - # output structures - cb = list(itertools.combinations(range(num_cell), 2)) - - for i in range(len(cb)): - arr1 = np.stack([x, mat[cb[i][0], cb[i][1], :]], axis=1) - arr2 = np.stack( - [ - x, - s_mean[cb[i][0], cb[i][1], :], - s_confint[cb[i][0], cb[i][1], :, 0], - s_confint[cb[i][0], cb[i][1], :, 1], - ], - axis=1, - ) - name = f"{str(cb[i][0])}-{str(cb[i][1])}" - info[name] = TimeSeriesData(arr1.T, file_name=name) - name = f"shuffle {str(cb[i][0])}-{str(cb[i][1])}" - info[name] = TimeSeriesData(arr2.T, file_name=name) + # output structures + cb = list(itertools.combinations(range(num_cell), 2)) + + for i in range(len(cb)): + arr1 = np.stack([x, mat[cb[i][0], cb[i][1], :]], axis=1) + arr2 = np.stack( + [ + x, + s_mean[cb[i][0], cb[i][1], :], + s_confint[cb[i][0], cb[i][1], :, 0], + s_confint[cb[i][0], cb[i][1], :, 1], + ], + axis=1, + ) + + name = f"{str(cb[i][0])}-{str(cb[i][1])}" + info[name] = TimeSeriesData(arr1.T, file_name=name) + name = f"shuffle {str(cb[i][0])}-{str(cb[i][1])}" + info[name] = TimeSeriesData(arr2.T, file_name=name) - return info + return info diff --git a/studio/app/optinist/wrappers/optinist/neural_population_analysis/granger.py b/studio/app/optinist/wrappers/optinist/neural_population_analysis/granger.py index c77597164..8afeea878 100644 --- a/studio/app/optinist/wrappers/optinist/neural_population_analysis/granger.py +++ b/studio/app/optinist/wrappers/optinist/neural_population_analysis/granger.py @@ -1,179 +1,221 @@ +from typing import Iterable, Optional, Union + +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import HeatMapData, ScatterData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import FluoData, IscellData from studio.app.optinist.wrappers.optinist.utils import standard_norm -def Granger( - neural_data: FluoData, - output_dir: str, - iscell: IscellData = None, - params: dict = None, - **kwargs, -) -> dict(): - # modules specific to function - # from sklearn.preprocessing import StandardScaler - import itertools - - import numpy as np - from statsmodels.tsa.stattools import adfuller, coint, grangercausalitytests - from tqdm import tqdm - - function_id = output_dir.split("/")[-1] - print("start granger:", function_id) - - neural_data = neural_data.data - - # data should be time x component matrix - if params["transpose"]: - X = neural_data.transpose() - else: - X = neural_data - - if iscell is not None: - iscell = iscell.data - ind = np.where(iscell > 0)[0] - X = X[:, ind] - - num_cell = X.shape[1] - comb = list(itertools.permutations(range(num_cell), 2)) # combinations with dup - num_comb = len(comb) - - # preprocessing - tX = standard_norm(X, params["standard_mean"], params["standard_std"]) - - # calculate dickey-fuller test - # augmented dickey-fuller test - # - if p val is large - # -> it cannot reject there is a unit root - # - small p-val means OK - # -> means this is not unit root process that it can apply Causality test - - adf = { - "adf_teststat": np.zeros([num_cell], dtype="float64"), - "adf_pvalue": np.zeros([num_cell], dtype="float64"), - "adf_usedlag": np.zeros([num_cell], dtype="int"), - "adf_nobs": np.zeros([num_cell], dtype="int"), - "adf_critical_values": np.zeros([num_cell, 3], dtype="float64"), - "adf_icbest": np.zeros([num_cell], dtype="float64"), - } - - if params["use_adfuller_test"]: - print("adfuller test ") - - for i in tqdm(range(num_cell)): - tp = adfuller(tX[:, i], **params["adfuller"]) - - adf["adf_teststat"][i] = tp[0] - adf["adf_pvalue"][i] = tp[1] - adf["adf_usedlag"][i] = tp[2] - adf["adf_nobs"][i] = tp[3] - adf["adf_critical_values"][i, :] = np.array( - [tp[4]["1%"], tp[4]["5%"], tp[4]["10%"]] - ) - adf["adf_icbest"][i] = tp[5] - - # test for cointegration - # augmented engle-granger two-step test - # Test for no-cointegration of a univariate equation - # if p val is small, the relation is cointegration - # -> check this if ADF pval is large - cit = { - "cit_count_t": np.zeros([num_comb], dtype="float64"), - "cit_pvalue": np.zeros([num_comb], dtype="int"), - "cit_crit_value": np.zeros([num_comb, 3], dtype="float64"), - } - - if params["use_coint_test"]: - print("cointegration test ") - - for i in tqdm(range(num_comb)): - tp = coint(X[:, comb[i][0]], X[:, comb[i][1]], **params["coint"]) - if not np.isnan(tp[0]): - cit["cit_count_t"][i] = tp[0] - if not np.isnan(tp[1]): - cit["cit_pvalue"][i] = tp[1] - cit["cit_crit_value"][i, :] = tp[2] - - # Granger causality - print("granger test ") - - if hasattr(params["Granger_maxlag"], "__iter__"): - num_lag = len(params["Granger_maxlag"]) - else: - num_lag = params["Granger_maxlag"] - - GC = { - "gc_combinations": comb, - "gc_ssr_ftest": np.zeros([num_comb, num_lag, 4], dtype="float64"), - "gc_ssr_chi2test": np.zeros([num_comb, num_lag, 3], dtype="float64"), - "gc_lrtest": np.zeros([num_comb, num_lag, 3], dtype="float64"), - "gc_params_ftest": np.zeros([num_comb, num_lag, 4], dtype="float64"), - "gc_OLS_restricted": [[0] * num_lag for i in range(num_comb)], - "gc_OLS_unrestricted": [[0] * num_lag for i in range(num_comb)], - "gc_OLS_restriction_matrix": [[0] * num_lag for i in range(num_comb)], - "Granger_fval_mat": [np.zeros([num_cell, num_cell]) for i in range(num_lag)], - } - - for i in tqdm(range(len(comb))): - # The Null hypothesis for grangercausalitytests is - # that the time series in the second column1, - # does NOT Granger cause the time series in the first column0 - # column 1 -> column 0 - tp = grangercausalitytests( - tX[:, [comb[i][0], comb[i][1]]], - params["Granger_maxlag"], - verbose=False, - addconst=params["Granger_addconst"], - ) +class GrangerAdfullerParams(BaseModel): + maxlag: Optional[int] = Field(None) + regression: str = Field("c") + autolag: str = Field("AIC") + store: bool = Field(False) + regresults: bool = Field(False) + - for j in range(len(tp)): # number of lag - GC["gc_ssr_ftest"][i, j, :] = tp[j + 1][0]["ssr_ftest"][ - 0:4 - ] # ssr based F test (F, pval, df_denom, df_num) - GC["gc_ssr_chi2test"][i, j, :] = tp[j + 1][0]["ssr_chi2test"][ - 0:3 - ] # ssr based chi2test (chi2, pval, df) - GC["gc_lrtest"][i, j, :] = tp[j + 1][0]["lrtest"][ - 0:3 - ] # likelihood ratio test (chi2, pval, df) - GC["gc_params_ftest"][i, j, :] = tp[j + 1][0]["params_ftest"][ - 0:4 - ] # parameter F test (F, pval, df_denom, df_num) - GC["gc_OLS_restricted"][i][j] = tp[j + 1][1][0] - GC["gc_OLS_unrestricted"][i][j] = tp[j + 1][1][1] - GC["gc_OLS_restriction_matrix"][i][j] = tp[j + 1][1][2] - - GC["Granger_fval_mat"][j][comb[i][0], comb[i][1]] = tp[j + 1][0][ - "ssr_ftest" - ][0] - - GC["Granger_fval_mat"] = np.array(GC["Granger_fval_mat"]) - - # main results for plot - info = {} - info["Granger_fval_mat_heatmap"] = HeatMapData( - GC["Granger_fval_mat"][0], file_name="gfm_heatmap" +class GrangerCointParams(BaseModel): + trend: str = Field("c") + method: str = Field("aeg") + maxlag: Optional[int] = Field(None) + autolag: str = Field("AIC") + + +class GrangerParams(BaseModel): + standard_mean: bool = Field(True) + standard_std: bool = Field(True) + transpose: bool = Field(True) + Granger_maxlag: Union[int, Iterable[int]] = Field( + 1, + description=( + "If an integer, computes the test for all lags up to maxlag. " + "If an iterable, computes the tests only for the lags in maxlag." + ), ) - info["Granger_fval_mat_scatter"] = ScatterData( - GC["Granger_fval_mat"][0], file_name="gfm" + Granger_addconst: bool = Field( + True, description="Include a constant in the model or not." ) + use_adfuller_test: bool = Field(True) + use_coint_test: bool = Field(True) + adfuller: GrangerAdfullerParams = Field(GrangerAdfullerParams()) + coint: GrangerCointParams = Field(GrangerCointParams()) + + +class Granger(AlgoTemplate): + def run( + self, + params: GrangerParams, + neural_data: FluoData, + iscell: IscellData = None, + ) -> dict(): + # modules specific to function + # from sklearn.preprocessing import StandardScaler + import itertools + + import numpy as np + from statsmodels.tsa.stattools import adfuller, coint, grangercausalitytests + from tqdm import tqdm + + print("start granger:", self.function_id) + + neural_data = neural_data.data + + # data should be time x component matrix + if params["transpose"]: + X = neural_data.transpose() + else: + X = neural_data + + if iscell is not None: + iscell = iscell.data + ind = np.where(iscell > 0)[0] + X = X[:, ind] + + num_cell = X.shape[1] + comb = list(itertools.permutations(range(num_cell), 2)) # combinations with dup + num_comb = len(comb) + + # preprocessing + tX = standard_norm(X, params["standard_mean"], params["standard_std"]) + + # calculate dickey-fuller test + # augmented dickey-fuller test + # - if p val is large + # -> it cannot reject there is a unit root + # - small p-val means OK + # -> means this is not unit root process that it can apply Causality test + + adf = { + "adf_teststat": np.zeros([num_cell], dtype="float64"), + "adf_pvalue": np.zeros([num_cell], dtype="float64"), + "adf_usedlag": np.zeros([num_cell], dtype="int"), + "adf_nobs": np.zeros([num_cell], dtype="int"), + "adf_critical_values": np.zeros([num_cell, 3], dtype="float64"), + "adf_icbest": np.zeros([num_cell], dtype="float64"), + } + + if params["use_adfuller_test"]: + print("adfuller test ") + + for i in tqdm(range(num_cell)): + tp = adfuller(tX[:, i], **params["adfuller"]) + + adf["adf_teststat"][i] = tp[0] + adf["adf_pvalue"][i] = tp[1] + adf["adf_usedlag"][i] = tp[2] + adf["adf_nobs"][i] = tp[3] + adf["adf_critical_values"][i, :] = np.array( + [tp[4]["1%"], tp[4]["5%"], tp[4]["10%"]] + ) + adf["adf_icbest"][i] = tp[5] + + # test for cointegration + # augmented engle-granger two-step test + # Test for no-cointegration of a univariate equation + # if p val is small, the relation is cointegration + # -> check this if ADF pval is large + cit = { + "cit_count_t": np.zeros([num_comb], dtype="float64"), + "cit_pvalue": np.zeros([num_comb], dtype="int"), + "cit_crit_value": np.zeros([num_comb, 3], dtype="float64"), + } + + if params["use_coint_test"]: + print("cointegration test ") + + for i in tqdm(range(num_comb)): + tp = coint(X[:, comb[i][0]], X[:, comb[i][1]], **params["coint"]) + if not np.isnan(tp[0]): + cit["cit_count_t"][i] = tp[0] + if not np.isnan(tp[1]): + cit["cit_pvalue"][i] = tp[1] + cit["cit_crit_value"][i, :] = tp[2] + + # Granger causality + print("granger test ") + + if hasattr(params["Granger_maxlag"], "__iter__"): + num_lag = len(params["Granger_maxlag"]) + else: + num_lag = params["Granger_maxlag"] + + GC = { + "gc_combinations": comb, + "gc_ssr_ftest": np.zeros([num_comb, num_lag, 4], dtype="float64"), + "gc_ssr_chi2test": np.zeros([num_comb, num_lag, 3], dtype="float64"), + "gc_lrtest": np.zeros([num_comb, num_lag, 3], dtype="float64"), + "gc_params_ftest": np.zeros([num_comb, num_lag, 4], dtype="float64"), + "gc_OLS_restricted": [[0] * num_lag for i in range(num_comb)], + "gc_OLS_unrestricted": [[0] * num_lag for i in range(num_comb)], + "gc_OLS_restriction_matrix": [[0] * num_lag for i in range(num_comb)], + "Granger_fval_mat": [ + np.zeros([num_cell, num_cell]) for i in range(num_lag) + ], + } + + for i in tqdm(range(len(comb))): + # The Null hypothesis for grangercausalitytests is + # that the time series in the second column1, + # does NOT Granger cause the time series in the first column0 + # column 1 -> column 0 + tp = grangercausalitytests( + tX[:, [comb[i][0], comb[i][1]]], + params["Granger_maxlag"], + verbose=False, + addconst=params["Granger_addconst"], + ) + + for j in range(len(tp)): # number of lag + GC["gc_ssr_ftest"][i, j, :] = tp[j + 1][0]["ssr_ftest"][ + 0:4 + ] # ssr based F test (F, pval, df_denom, df_num) + GC["gc_ssr_chi2test"][i, j, :] = tp[j + 1][0]["ssr_chi2test"][ + 0:3 + ] # ssr based chi2test (chi2, pval, df) + GC["gc_lrtest"][i, j, :] = tp[j + 1][0]["lrtest"][ + 0:3 + ] # likelihood ratio test (chi2, pval, df) + GC["gc_params_ftest"][i, j, :] = tp[j + 1][0]["params_ftest"][ + 0:4 + ] # parameter F test (F, pval, df_denom, df_num) + GC["gc_OLS_restricted"][i][j] = tp[j + 1][1][0] + GC["gc_OLS_unrestricted"][i][j] = tp[j + 1][1][1] + GC["gc_OLS_restriction_matrix"][i][j] = tp[j + 1][1][2] + + GC["Granger_fval_mat"][j][comb[i][0], comb[i][1]] = tp[j + 1][0][ + "ssr_ftest" + ][0] + + GC["Granger_fval_mat"] = np.array(GC["Granger_fval_mat"]) + + # main results for plot + info = {} + info["Granger_fval_mat_heatmap"] = HeatMapData( + GC["Granger_fval_mat"][0], file_name="gfm_heatmap" + ) + info["Granger_fval_mat_scatter"] = ScatterData( + GC["Granger_fval_mat"][0], file_name="gfm" + ) - # NWB追加 - nwbfile = {} - nwbfile[NWBDATASET.POSTPROCESS] = { - function_id: { - "Granger_fval_mat": GC["Granger_fval_mat"][0], - "gc_combinations": GC["gc_combinations"], - "gc_ssr_ftest": GC["gc_ssr_ftest"], - "gc_ssr_chi2test": GC["gc_ssr_chi2test"], - "gc_lrtest": GC["gc_lrtest"], - "gc_params_ftest": GC["gc_params_ftest"], - "cit_pvalue": cit["cit_pvalue"], - "adf_pvalue": adf["adf_pvalue"], + # NWB追加 + nwbfile = {} + nwbfile[NWBDATASET.POSTPROCESS] = { + self.function_id: { + "Granger_fval_mat": GC["Granger_fval_mat"][0], + "gc_combinations": GC["gc_combinations"], + "gc_ssr_ftest": GC["gc_ssr_ftest"], + "gc_ssr_chi2test": GC["gc_ssr_chi2test"], + "gc_lrtest": GC["gc_lrtest"], + "gc_params_ftest": GC["gc_params_ftest"], + "cit_pvalue": cit["cit_pvalue"], + "adf_pvalue": adf["adf_pvalue"], + } } - } - info["nwbfile"] = nwbfile + info["nwbfile"] = nwbfile - return info + return info diff --git a/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/correlation.yaml b/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/correlation.yaml deleted file mode 100644 index ca5f01de0..000000000 --- a/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/correlation.yaml +++ /dev/null @@ -1 +0,0 @@ -transpose: True diff --git a/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/cross_correlation.yaml b/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/cross_correlation.yaml deleted file mode 100644 index 458700cc5..000000000 --- a/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/cross_correlation.yaml +++ /dev/null @@ -1,13 +0,0 @@ -transpose: False - -method: 'direct' - -# lags: int number of frames to show (+-frames) -lags: 100 - -# shuffle_sample_number: int -# number of cross correlation caliculation for randomization -shuffle_sample_number: 100 - -# shuffle_confidence_interval: float -shuffle_confidence_interval: 0.95 diff --git a/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/granger.yaml b/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/granger.yaml deleted file mode 100644 index 3f66a68a4..000000000 --- a/studio/app/optinist/wrappers/optinist/neural_population_analysis/params/granger.yaml +++ /dev/null @@ -1,27 +0,0 @@ -standard_mean: True -standard_std: True - -transpose: True - -# Granger_maxlag(int, Iterable[int]): If an integer, computes the test for all lags up to maxlag. If an iterable, computes the tests only for the lags in maxlag. -# no default value -Granger_maxlag: 1 - -#Include a constant in the model or not. -Granger_addconst: True - -use_adfuller_test: True -use_coint_test: True - -adfuller: - maxlag: - regression: 'c' - autolag: 'AIC' - store: False - regresults: False - -coint: - trend: 'c' - method: 'aeg' - maxlag: - autolag: 'AIC' diff --git a/studio/app/optinist/wrappers/suite2p/__init__.py b/studio/app/optinist/wrappers/suite2p/__init__.py index 07243b7b1..533d833bf 100644 --- a/studio/app/optinist/wrappers/suite2p/__init__.py +++ b/studio/app/optinist/wrappers/suite2p/__init__.py @@ -1,24 +1,24 @@ -from studio.app.optinist.wrappers.suite2p.file_convert import suite2p_file_convert -from studio.app.optinist.wrappers.suite2p.registration import suite2p_registration -from studio.app.optinist.wrappers.suite2p.roi import suite2p_roi -from studio.app.optinist.wrappers.suite2p.spike_deconv import suite2p_spike_deconv +from studio.app.optinist.wrappers.suite2p.file_convert import Suite2pFileConvert +from studio.app.optinist.wrappers.suite2p.registration import Suite2pRegistration +from studio.app.optinist.wrappers.suite2p.roi import Suite2pRoi +from studio.app.optinist.wrappers.suite2p.spike_deconv import Suite2pSpikeDeconv suite2p_wrapper_dict = { "suite2p": { "suite2p_file_convert": { - "function": suite2p_file_convert, + "function": Suite2pFileConvert, "conda_name": "suite2p", }, "suite2p_registration": { - "function": suite2p_registration, + "function": Suite2pRegistration, "conda_name": "suite2p", }, "suite2p_roi": { - "function": suite2p_roi, + "function": Suite2pRoi, "conda_name": "suite2p", }, "suite2p_spike_deconv": { - "function": suite2p_spike_deconv, + "function": Suite2pSpikeDeconv, "conda_name": "suite2p", }, } diff --git a/studio/app/optinist/wrappers/suite2p/file_convert.py b/studio/app/optinist/wrappers/suite2p/file_convert.py index ec183ea5f..ec8baebb8 100644 --- a/studio/app/optinist/wrappers/suite2p/file_convert.py +++ b/studio/app/optinist/wrappers/suite2p/file_convert.py @@ -1,5 +1,9 @@ import os +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.core.utils.filepath_creater import ( create_directory, join_filepath, @@ -8,43 +12,51 @@ from studio.app.optinist.dataclass import Suite2pData -def suite2p_file_convert( - image: ImageData, output_dir: str, params: dict = None, **kwargs -) -> dict(ops=Suite2pData): - from suite2p import default_ops, io +class Suite2pFileConvertParams(BaseModel): + nplanes: int = Field(1) + nchannels: int = Field(1) + force_sktiff: bool = Field(False) + batch_size: int = Field(500) + do_registration: int = Field(1) + + +class Suite2pFileConvert(AlgoTemplate): + def run( + self, params: Suite2pFileConvertParams, image: ImageData + ) -> dict(ops=Suite2pData): + from suite2p import default_ops, io - function_id = output_dir.split("/")[-1] - print("start suite2p_file_convert:", function_id) + print("start suite2p_file_convert:", self.function_id) - data_path_list = [] - data_name_list = [] - for file_path in image.path: - data_path_list.append(os.path.dirname(file_path)) - data_name_list.append(os.path.basename(file_path)) + data_path_list = [] + data_name_list = [] + for file_path in image.path: + data_path_list.append(os.path.dirname(file_path)) + data_name_list.append(os.path.basename(file_path)) - print(data_path_list) - print(data_name_list) - # data pathと保存pathを指定 - db = { - "data_path": data_path_list, - "tiff_list": data_name_list, - "save_path0": output_dir, - "save_folder": "suite2p", - } + print(data_path_list) + print(data_name_list) + # data pathと保存pathを指定 + db = { + "data_path": data_path_list, + "tiff_list": data_name_list, + "save_path0": self.output_dir, + "save_folder": "suite2p", + } - ops = {**default_ops(), **params, **db} + ops = {**default_ops(), **params, **db} - # save folderを指定 - create_directory(join_filepath([ops["save_path0"], ops["save_folder"]])) + # save folderを指定 + create_directory(join_filepath([ops["save_path0"], ops["save_folder"]])) - # save ops.npy(parameter) and data.bin - ops = io.tiff_to_binary(ops.copy()) + # save ops.npy(parameter) and data.bin + ops = io.tiff_to_binary(ops.copy()) - info = { - "meanImg": ImageData( - ops["meanImg"], output_dir=output_dir, file_name="meanImg" - ), - "ops": Suite2pData(ops, file_name="ops"), - } + info = { + "meanImg": ImageData( + ops["meanImg"], output_dir=self.output_dir, file_name="meanImg" + ), + "ops": Suite2pData(ops, file_name="ops"), + } - return info + return info diff --git a/studio/app/optinist/wrappers/suite2p/params/suite2p_file_convert.yaml b/studio/app/optinist/wrappers/suite2p/params/suite2p_file_convert.yaml deleted file mode 100644 index 3766c078f..000000000 --- a/studio/app/optinist/wrappers/suite2p/params/suite2p_file_convert.yaml +++ /dev/null @@ -1,5 +0,0 @@ -nplanes: 1 -nchannels: 1 -force_sktiff: False -batch_size: 500 -do_registration: 1 diff --git a/studio/app/optinist/wrappers/suite2p/params/suite2p_registration.yaml b/studio/app/optinist/wrappers/suite2p/params/suite2p_registration.yaml deleted file mode 100644 index afce2ac66..000000000 --- a/studio/app/optinist/wrappers/suite2p/params/suite2p_registration.yaml +++ /dev/null @@ -1,24 +0,0 @@ -frames_include: -1 -keep_movie_raw: False -do_bidiphase: False - -smooth_sigma: 1.15 -smooth_sigma_time: 0 -bidiphase: 0 -maxregshift: 0.1 -maxregshiftNR: 5 -nonrigid: True -block_size: [128, 128] -snr_thresh: 1.2 -functional_chan : 1 -align_by_chan : 1 -reg_tif: False -th_badframes: 1.0 -diameter: 0 - -# 1P setting -1Preg: False -spatial_hp_reg: 42 -pre_smooth: 0 -spatial_taper: 40 -bidi_corrected: False diff --git a/studio/app/optinist/wrappers/suite2p/params/suite2p_roi.yaml b/studio/app/optinist/wrappers/suite2p/params/suite2p_roi.yaml deleted file mode 100644 index 8d4c731c7..000000000 --- a/studio/app/optinist/wrappers/suite2p/params/suite2p_roi.yaml +++ /dev/null @@ -1,28 +0,0 @@ -# main settings -tau: 1.0 # this is the main parameter for deconvolution - -# classification parameters -soma_crop: True # crop dendrites for cell classification stats like compactness - -# cell detection settings -high_pass: 100 # running mean subtraction with window of size 'high_pass(use low values for 1P) -sparse_mode: True # whether or not to run sparse_mode -max_overlap: 0.75 # cells with more overlap than this get removed during triage before refinement -nbinned: 5000 # max number of binned frames for cell detection -spatial_scale: 0 # 0: multi-scale; 1: 6 pixels 2: 12 pixels 3: 24 pixels 4: 48 pixels -threshold_scaling: 1.0 # adjust the automatically determined threshold by this scalar multiplier -max_iterations: 20 # maximum number of iterations to do cell detection - -# 1P settings -spatial_hp_detect: 25 # window for spatial high-pass filtering for neuropil subtraction before detection - -# output settings -preclassify: 0. # apply classifier before signal extraction with probability 0.3 - -# ROI extraction parameters -allow_overlap: False # pixels that are overlapping are thrown out (False) or added to both ROIs (True) -inner_neuropil_radius: 2 # number of pixels to keep between ROI and neuropil donut -min_neuropil_pixels: 350 # minimum number of pixels in the neuropil - -# deconvolution settings -neucoeff: .7 # neuropil coefficient diff --git a/studio/app/optinist/wrappers/suite2p/params/suite2p_spike_deconv.yaml b/studio/app/optinist/wrappers/suite2p/params/suite2p_spike_deconv.yaml deleted file mode 100644 index 37c1973d9..000000000 --- a/studio/app/optinist/wrappers/suite2p/params/suite2p_spike_deconv.yaml +++ /dev/null @@ -1,6 +0,0 @@ -# deconvolution settings -baseline: 'maximin' # baselining mode (can also choose 'prctile') -win_baseline: 60. # window for maximin -sig_baseline: 10. # smoothing constant for gaussian filter -prctile_baseline: 8. # optional (whether to use a percentile baseline) -neucoeff: .7 # neuropil coefficient diff --git a/studio/app/optinist/wrappers/suite2p/registration.py b/studio/app/optinist/wrappers/suite2p/registration.py index 9c5b64ba4..b435967b1 100644 --- a/studio/app/optinist/wrappers/suite2p/registration.py +++ b/studio/app/optinist/wrappers/suite2p/registration.py @@ -1,38 +1,71 @@ +from typing import List + +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import ImageData from studio.app.optinist.dataclass import Suite2pData -def suite2p_registration( - ops: Suite2pData, output_dir: str, params: dict = None, **kwargs -) -> dict(ops=Suite2pData): - from suite2p import default_ops, registration +class Suite2pRegistrationParams(BaseModel): + frames_include: int = Field(-1) + keep_movie_raw: bool = Field(False) + do_bidiphase: bool = Field(False) + + smooth_sigma: float = Field(1.15) + smooth_sigma_time: int = Field(0) + bidiphase: int = Field(0) + maxregshift: float = Field(0.1) + maxregshiftNR: int = Field(5) + nonrigid: bool = Field(True) + block_size: List[int] = Field([128, 128]) + snr_thresh: float = Field(1.2) + functional_chan: int = Field(1) + align_by_chan: int = Field(1) + reg_tif: bool = Field(False) + th_badframes: float = Field(1.0) + diameter: int = Field(0) + one_preg: bool = Field(False, alias="1Preg") + spatial_hp_reg: int = Field(42) + pre_smooth: int = Field(0) + spatial_taper: int = Field(40) + bidi_corrected: bool = Field(False) + + +class Suite2pRegistration(AlgoTemplate): + def run( + self, params: Suite2pRegistrationParams, ops: Suite2pData + ) -> dict(ops=Suite2pData): + from suite2p import default_ops, registration - function_id = output_dir.split("/")[-1] - print("start suite2p registration:", function_id) + print("start suite2p registration:", self.function_id) - ops = ops.data - refImg = ops["meanImg"] - print("start suite2_registration") + ops = ops.data + refImg = ops["meanImg"] + print("start suite2_registration") - # REGISTRATION - if len(refImg.shape) == 3: - refImg = refImg[0] + # REGISTRATION + if len(refImg.shape) == 3: + refImg = refImg[0] - ops = {**default_ops(), **ops, **params} + ops = {**default_ops(), **ops, **params} - # register binary - ops = registration.register_binary(ops, refImg=refImg) + # register binary + ops = registration.register_binary(ops, refImg=refImg) - # compute metrics for registration - if ops.get("do_regmetrics", True) and ops["nframes"] >= 1500: - ops = registration.get_pc_metrics(ops) + # compute metrics for registration + if ops.get("do_regmetrics", True) and ops["nframes"] >= 1500: + ops = registration.get_pc_metrics(ops) - info = { - "refImg": ImageData(ops["refImg"], output_dir=output_dir, file_name="refImg"), - "meanImgE": ImageData( - ops["meanImgE"], output_dir=output_dir, file_name="meanImgE" - ), - "ops": Suite2pData(ops, file_name="ops"), - } + info = { + "refImg": ImageData( + ops["refImg"], output_dir=self.output_dir, file_name="refImg" + ), + "meanImgE": ImageData( + ops["meanImgE"], output_dir=self.output_dir, file_name="meanImgE" + ), + "ops": Suite2pData(ops, file_name="ops"), + } - return info + return info diff --git a/studio/app/optinist/wrappers/suite2p/roi.py b/studio/app/optinist/wrappers/suite2p/roi.py index 3c077ca16..ea8c91f37 100644 --- a/studio/app/optinist/wrappers/suite2p/roi.py +++ b/studio/app/optinist/wrappers/suite2p/roi.py @@ -1,130 +1,219 @@ +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.common.dataclass import ImageData from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import FluoData, IscellData, RoiData, Suite2pData -def suite2p_roi( - ops: Suite2pData, output_dir: str, params: dict = None, **kwargs -) -> dict(ops=Suite2pData, fluorescence=FluoData, iscell=IscellData): - import numpy as np - from suite2p import ROI, classification, default_ops, detection, extraction - - function_id = output_dir.split("/")[-1] - print("start suite2p_roi:", function_id) - - nwbfile = kwargs.get("nwbfile", {}) - fs = nwbfile.get("imaging_plane", {}).get("imaging_rate", 30) - - ops = ops.data - ops = {**default_ops(), **ops, **params, "fs": fs} - - # ROI detection - ops_classfile = ops.get("classifier_path") - builtin_classfile = classification.builtin_classfile - user_classfile = classification.user_classfile - if ops_classfile: - print(f"NOTE: applying classifier {str(ops_classfile)}") - classfile = ops_classfile - elif ops["use_builtin_classifier"] or not user_classfile.is_file(): - print(f"NOTE: Applying builtin classifier at {str(builtin_classfile)}") - classfile = builtin_classfile - else: - print(f"NOTE: applying default {str(user_classfile)}") - classfile = user_classfile - - ops, stat = detection.detect(ops=ops, classfile=classfile) - - # ROI EXTRACTION - ops, stat, F, Fneu, _, _ = extraction.create_masks_and_extract(ops, stat) - stat = stat.tolist() - - # ROI CLASSIFICATION - iscell = classification.classify(stat=stat, classfile=classfile) - iscell = iscell[:, 0].astype(bool) - - arrays = [] - for i, s in enumerate(stat): - array = ROI( - ypix=s["ypix"], xpix=s["xpix"], lam=s["lam"], med=s["med"], do_crop=False - ).to_array(Ly=ops["Ly"], Lx=ops["Lx"]) - array *= i + 1 - arrays.append(array) - - im = np.stack(arrays) - im[im == 0] = np.nan - im -= 1 - - # roiを追加 - roi_list = [] - for i in range(len(stat)): - kargs = {} - kargs["pixel_mask"] = np.array( - [stat[i]["ypix"], stat[i]["xpix"], stat[i]["lam"]] - ).T - roi_list.append(kargs) - - # NWBを追加 - nwbfile = {} - - nwbfile[NWBDATASET.ROI] = {function_id: roi_list} - - # iscellを追加 - nwbfile[NWBDATASET.COLUMN] = { - function_id: { - "name": "iscell", - "description": "two columns - iscell & probcell", - "data": iscell, - } - } - - # Fluorenceを追加 - nwbfile[NWBDATASET.FLUORESCENCE] = { - function_id: { - "Fluorescence": { - "table_name": "Fluorescence", - "region": list(range(len(F))), - "name": "Fluorescence", - "data": F, - "unit": "lumens", - "rate": ops["fs"], - }, - "Neuropil": { - "table_name": "Neuropil", - "region": list(range(len(Fneu))), - "name": "Neuropil", - "data": Fneu, - "unit": "lumens", - "rate": ops["fs"], - }, - } - } - - ops["stat"] = stat - ops["F"] = F - ops["Fneu"] = Fneu - ops["iscell"] = iscell - ops["im"] = im - - info = { - "ops": Suite2pData(ops), - "max_proj": ImageData( - ops["max_proj"], output_dir=output_dir, file_name="max_proj" +class Suite2pRoiParams(BaseModel): + tau: float = Field(1.0, description="this is the main parameter for deconvolution") + + # classification parameters + soma_crop: bool = Field( + True, + description="crop dendrites for cell classification stats like compactness", + ) + + # cell detection settings + high_pass: int = Field( + 100, + description=( + "running mean subtraction with window of size " + "'high_pass(use low values for 1P" + ), + ) + sparse_mode: bool = Field(True, description="whether or not to run sparse_mode") + max_overlap: float = Field( + 0.75, + description=( + "cells with more overlap than this " + "get removed during triage before refinement" ), - "Vcorr": ImageData(ops["Vcorr"], output_dir=output_dir, file_name="Vcorr"), - "fluorescence": FluoData(F, file_name="fluorescence"), - "iscell": IscellData(iscell, file_name="iscell"), - "all_roi": RoiData( - np.nanmax(im, axis=0), output_dir=output_dir, file_name="all_roi" + ) + nbinned: int = Field( + 5000, description=" max number of binned frames for cell detection" + ) + spatial_scale: int = Field( + 0, + description=( + "0: multi-scale; 1: 6 pixels " "2: 12 pixels 3: 24 pixels 4: 48 pixels" ), - "non_cell_roi": RoiData( - np.nanmax(im[~iscell], axis=0), - output_dir=output_dir, - file_name="noncell_roi", + ) + threshold_scaling: float = Field( + 1.0, + description=( + "adjust the automatically determined threshold " "by this scalar multiplier" ), - "cell_roi": RoiData( - np.nanmax(im[iscell], axis=0), output_dir=output_dir, file_name="cell_roi" + ) + max_iterations: int = Field( + 20, description="maximum number of iterations to do cell detection" + ) + + # 1P settings + spatial_hp_detect: int = Field( + 25, + description=( + "window for spatial high-pass filtering for neuropil subtraction " + "before detection" ), - "nwbfile": nwbfile, - } + ) + + # output settings + preclassify: float = Field( + 0.0, + description="apply classifier before signal extraction with probability 0.3", + ) + + # ROI extraction parameters + allow_overlap: bool = Field( + False, + description=( + "pixels that are overlapping are thrown out (False) " + "or added to both ROIs (True)" + ), + ) + inner_neuropil_radius: int = Field( + 2, description="number of pixels to keep between ROI and neuropil donut" + ) + min_neuropil_pixels: int = Field( + 350, description="minimum number of pixels in the neuropil" + ) + + # deconvolution settings + neucoeff: float = Field(0.7, description="neuropil coefficient") + + +class Suite2pRoi(AlgoTemplate): + def run( + self, params: Suite2pRoiParams, ops: Suite2pData + ) -> dict(ops=Suite2pData, fluorescence=FluoData, iscell=IscellData): + import numpy as np + from suite2p import ROI, classification, default_ops, detection, extraction + + print("start suite2p_roi:", self.function_id) + + fs = self.nwb_params.get("imaging_plane", {}).get("imaging_rate", 30) + + ops = ops.data + ops = {**default_ops(), **ops, **params, "fs": fs} + + # ROI detection + ops_classfile = ops.get("classifier_path") + builtin_classfile = classification.builtin_classfile + user_classfile = classification.user_classfile + if ops_classfile: + print(f"NOTE: applying classifier {str(ops_classfile)}") + classfile = ops_classfile + elif ops["use_builtin_classifier"] or not user_classfile.is_file(): + print(f"NOTE: Applying builtin classifier at {str(builtin_classfile)}") + classfile = builtin_classfile + else: + print(f"NOTE: applying default {str(user_classfile)}") + classfile = user_classfile + + ops, stat = detection.detect(ops=ops, classfile=classfile) + + # ROI EXTRACTION + ops, stat, F, Fneu, _, _ = extraction.create_masks_and_extract(ops, stat) + stat = stat.tolist() + + # ROI CLASSIFICATION + iscell = classification.classify(stat=stat, classfile=classfile) + iscell = iscell[:, 0].astype(bool) + + arrays = [] + for i, s in enumerate(stat): + array = ROI( + ypix=s["ypix"], + xpix=s["xpix"], + lam=s["lam"], + med=s["med"], + do_crop=False, + ).to_array(Ly=ops["Ly"], Lx=ops["Lx"]) + array *= i + 1 + arrays.append(array) + + im = np.stack(arrays) + im[im == 0] = np.nan + im -= 1 + + # roiを追加 + roi_list = [] + for i in range(len(stat)): + kargs = {} + kargs["pixel_mask"] = np.array( + [stat[i]["ypix"], stat[i]["xpix"], stat[i]["lam"]] + ).T + roi_list.append(kargs) + + # NWBを追加 + nwbfile = {} + + nwbfile[NWBDATASET.ROI] = {self.function_id: roi_list} + + # iscellを追加 + nwbfile[NWBDATASET.COLUMN] = { + self.function_id: { + "name": "iscell", + "description": "two columns - iscell & probcell", + "data": iscell, + } + } + + # Fluorenceを追加 + nwbfile[NWBDATASET.FLUORESCENCE] = { + self.function_id: { + "Fluorescence": { + "table_name": "Fluorescence", + "region": list(range(len(F))), + "name": "Fluorescence", + "data": F, + "unit": "lumens", + "rate": ops["fs"], + }, + "Neuropil": { + "table_name": "Neuropil", + "region": list(range(len(Fneu))), + "name": "Neuropil", + "data": Fneu, + "unit": "lumens", + "rate": ops["fs"], + }, + } + } + + ops["stat"] = stat + ops["F"] = F + ops["Fneu"] = Fneu + ops["iscell"] = iscell + ops["im"] = im + + info = { + "ops": Suite2pData(ops), + "max_proj": ImageData( + ops["max_proj"], output_dir=self.output_dir, file_name="max_proj" + ), + "Vcorr": ImageData( + ops["Vcorr"], output_dir=self.output_dir, file_name="Vcorr" + ), + "fluorescence": FluoData(F, file_name="fluorescence"), + "iscell": IscellData(iscell, file_name="iscell"), + "all_roi": RoiData( + np.nanmax(im, axis=0), output_dir=self.output_dir, file_name="all_roi" + ), + "non_cell_roi": RoiData( + np.nanmax(im[~iscell], axis=0), + output_dir=self.output_dir, + file_name="noncell_roi", + ), + "cell_roi": RoiData( + np.nanmax(im[iscell], axis=0), + output_dir=self.output_dir, + file_name="cell_roi", + ), + "nwbfile": nwbfile, + } - return info + return info diff --git a/studio/app/optinist/wrappers/suite2p/spike_deconv.py b/studio/app/optinist/wrappers/suite2p/spike_deconv.py index ad17592ee..07c248128 100644 --- a/studio/app/optinist/wrappers/suite2p/spike_deconv.py +++ b/studio/app/optinist/wrappers/suite2p/spike_deconv.py @@ -1,68 +1,86 @@ +from pydantic import BaseModel +from pydantic.dataclasses import Field + +from studio.app.common.core.algo import AlgoTemplate from studio.app.optinist.core.nwb.nwb import NWBDATASET from studio.app.optinist.dataclass import FluoData, Suite2pData -def suite2p_spike_deconv( - ops: Suite2pData, output_dir: str, params: dict = None, **kwargs -) -> dict(ops=Suite2pData, spks=FluoData): - import numpy as np - from suite2p import default_ops, extraction +class Suite2pSpikeDeconvParams(BaseModel): + baseline: str = Field( + "maximin", description="baselining mode (can also choose 'prctile')" + ) + win_baseline: float = Field(60.0, description="window for maximin") + sig_baseline: float = Field( + 10.0, description="smoothing constant for gaussian filter" + ) + prctile_baseline: float = Field( + 8.0, description="optional (whether to use a percentile baseline)" + ) + neucoeff: float = Field(0.7, description="neuropil coefficient") - function_id = output_dir.split("/")[-1] - print("start suite2_spike_deconv:", function_id) - ops = ops.data +class Suite2pSpikeDeconv(AlgoTemplate): + def run( + self, params: Suite2pSpikeDeconvParams, ops: Suite2pData + ) -> dict(ops=Suite2pData, spks=FluoData): + import numpy as np + from suite2p import default_ops, extraction - ops = {**default_ops(), **ops, **params} + print("start suite2_spike_deconv:", self.function_id) - dF = ops["F"].copy() - ops["neucoeff"] * ops["Fneu"] - dF = extraction.preprocess( - F=dF, - baseline=ops["baseline"], - win_baseline=ops["win_baseline"], - sig_baseline=ops["sig_baseline"], - fs=ops["fs"], - prctile_baseline=ops["prctile_baseline"], - ) - spks = extraction.oasis( - F=dF, batch_size=ops["batch_size"], tau=ops["tau"], fs=ops["fs"] - ) + ops = ops.data + + ops = {**default_ops(), **ops, **params} - ops["spks"] = spks - - # NWBを追加 - nwbfile = {} - - # roiを追加 - stat = ops["stat"] - roi_list = [] - for i in range(len(stat)): - kargs = {} - kargs["pixel_mask"] = np.array( - [stat[i]["ypix"], stat[i]["xpix"], stat[i]["lam"]] - ).T - roi_list.append(kargs) - - nwbfile[NWBDATASET.ROI] = {function_id: {"roi_list": roi_list}} - - # Fluorenceを追加 - nwbfile[NWBDATASET.FLUORESCENCE] = { - function_id: { - "Deconvolved": { - "table_name": "Deconvolved", - "region": list(range(len(spks))), - "name": function_id + "_Deconvolved", - "data": spks, - "unit": "lumens", - "rate": ops["fs"], + dF = ops["F"].copy() - ops["neucoeff"] * ops["Fneu"] + dF = extraction.preprocess( + F=dF, + baseline=ops["baseline"], + win_baseline=ops["win_baseline"], + sig_baseline=ops["sig_baseline"], + fs=ops["fs"], + prctile_baseline=ops["prctile_baseline"], + ) + spks = extraction.oasis( + F=dF, batch_size=ops["batch_size"], tau=ops["tau"], fs=ops["fs"] + ) + + ops["spks"] = spks + + # NWBを追加 + nwbfile = {} + + # roiを追加 + stat = ops["stat"] + roi_list = [] + for i in range(len(stat)): + kargs = {} + kargs["pixel_mask"] = np.array( + [stat[i]["ypix"], stat[i]["xpix"], stat[i]["lam"]] + ).T + roi_list.append(kargs) + + nwbfile[NWBDATASET.ROI] = {self.function_id: {"roi_list": roi_list}} + + # Fluorenceを追加 + nwbfile[NWBDATASET.FLUORESCENCE] = { + self.function_id: { + "Deconvolved": { + "table_name": "Deconvolved", + "region": list(range(len(spks))), + "name": self.function_id + "_Deconvolved", + "data": spks, + "unit": "lumens", + "rate": ops["fs"], + } } } - } - info = { - "ops": Suite2pData(ops), - "spks": FluoData(spks, file_name="spks"), - "nwbfile": nwbfile, - } + info = { + "ops": Suite2pData(ops), + "spks": FluoData(spks, file_name="spks"), + "nwbfile": nwbfile, + } - return info + return info diff --git a/studio/tests/app/common/core/snakemake/test_snakemake_setfile.py b/studio/tests/app/common/core/snakemake/test_snakemake_setfile.py index cfc5dec3a..9164aa7c1 100644 --- a/studio/tests/app/common/core/snakemake/test_snakemake_setfile.py +++ b/studio/tests/app/common/core/snakemake/test_snakemake_setfile.py @@ -8,7 +8,7 @@ id="input_0", type="type", data=NodeData( - label="label", + label="caiman_mc", param={}, path="path", type="type", diff --git a/studio/tests/app/common/core/utils/test_config_handler.py b/studio/tests/app/common/core/utils/test_config_handler.py index 35d481dd9..d2b9dad2f 100644 --- a/studio/tests/app/common/core/utils/test_config_handler.py +++ b/studio/tests/app/common/core/utils/test_config_handler.py @@ -1,30 +1,13 @@ import os -from studio.app.common.core.utils.config_handler import ConfigReader, ConfigWriter +from studio.app.common.core.utils.config_handler import ConfigWriter from studio.app.common.core.utils.filepath_creater import join_filepath -from studio.app.common.core.utils.filepath_finder import find_param_filepath from studio.app.dir_path import DIRPATH dirpath = f"{DIRPATH.DATA_DIR}/output" filename = "test.yaml" -def test_config_reader(): - filename = "eta" - filepath = find_param_filepath(filename) - config = ConfigReader.read(filepath) - - assert isinstance(config, dict) - assert len(config) > 0 - - filename = "not_exist_config" - filepath = find_param_filepath(filename) - config = ConfigReader.read(filepath) - - assert isinstance(config, dict) - assert len(config) == 0 - - def test_config_writer(): filepath = join_filepath([dirpath, filename]) diff --git a/studio/tests/app/common/routers/test_params.py b/studio/tests/app/common/routers/test_params.py index 26c5484a5..1702aed6b 100644 --- a/studio/tests/app/common/routers/test_params.py +++ b/studio/tests/app/common/routers/test_params.py @@ -5,11 +5,13 @@ def test_params(client): assert response.status_code == 200 assert isinstance(data, dict) - assert isinstance(data["border_nan"], str) - assert data["border_nan"] == "copy" + border_nan = data["border_nan"]["value"] + assert isinstance(border_nan, str) + assert border_nan == "copy" - assert isinstance(data["use_cuda"], bool) - assert data["use_cuda"] is False + use_cuda = data["use_cuda"]["value"] + assert isinstance(use_cuda, bool) + assert use_cuda is False def test_snakemake_params(client): @@ -19,8 +21,10 @@ def test_snakemake_params(client): assert response.status_code == 200 assert isinstance(data, dict) - assert isinstance(data["cores"], int) - assert data["cores"] == 2 + cores = data["cores"]["value"] + assert isinstance(cores, int) + assert cores == 2 - assert isinstance(data["use_conda"], bool) - assert data["use_conda"] is True + use_conda = data["use_conda"]["value"] + assert isinstance(use_conda, bool) + assert use_conda is True diff --git a/studio/tests/app/optinist/routers/test_nwb.py b/studio/tests/app/optinist/routers/test_nwb.py index db5380e20..1844d65c2 100644 --- a/studio/tests/app/optinist/routers/test_nwb.py +++ b/studio/tests/app/optinist/routers/test_nwb.py @@ -12,8 +12,9 @@ def test_nwb_params(client): assert response.status_code == 200 assert isinstance(data, dict) - assert isinstance(data["session_description"], str) - assert data["session_description"] == "optinist" + session_description = data["session_description"]["value"] + assert isinstance(session_description, str) + assert session_description == "optinist" def test_download_nwb(client):